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INTRODUCTION 

Large Liquid Rocket Engine Transient Performance Simulation. Contract 
NAS8-36994. is a program to design and develop a simulation system lor rocket 
engine transient performance models. The system was designed and 
conceptually demonstrated during Phase I of the contract. The current Phase lid 
the contract will enhance the system, develop a detailed simulation of the 
Technology Test Bed Engine (TTBE). and deliver the model and system isoWw* 

,o the National Aeronautics and Space Administration (NASA) George C_ 
Marshall Space Flight Center (MSFC) .The program is sponsored by NASA/MSFC 
“ W.A. Adams acting as project manager. The Pratt & Whitney program 

manager is Mr. J.R. Mason. 

The objective of the program is to apply state-of-the-art modeling 

simulation technology to the rocket engine models. This ap P ll0a 

to design develop, and produce all liquid rocket transient simulation system 

compatible with MSFC control analysis software. si ™' a, %' a ^. S R 
computers. The simulation system will be referred to as the ROCE S ( 

Engine Transient Simulation) System. 

The report covers the third six months of ROCETS development. 
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Phase 1 ot the ROCETS program consists of seven technical tasks plus provision 
for reports and review; these tasks are. 


1 . Architecture 

2. System Requirements 

3. Component and Submodel Requirements 

4. Submodel Implementation 


Component Implementation 
Submodel Testing and Verification 
Subsystem Testing and Verification 


Duri „ g Phase the above tasKs were combed and a ^ 
Phase I completion review were conducted at MS • 
of two technical tasks plus reports and reviews; these tas 


TTBE Model Data Generation 
System Testing and Verification 


During this period specific ooding of the system processors W* , and the 

engineering representations o, Phase 

-lure cen,^ on ^ 
variable common, GLOBVAR, were necessary to increase process 


The engineering modules completed during Phase II are listed below 


1 . INJTOO - Main Injector 

2. MCHBOO - Main Chamber 

3. NOZLOO - Nozzle Thrust Calculations 

4. PBRNOO - Preburner 

5 PIPE02 - Compressible flow without inertia 

6 PUMPOO - Polytropic Pump 

7 ROTR00 - Rotor Torque Balance/Speed Derivative 

8. TURBOO - Turbine 


Detailed documentation of these modules is in the appendix. addition to the 

submodules will be included in updates to the System Design sp 
(SDS). 
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Specific coding was begun on the system configuration processor. 
necessary for multiple module operation have been completed but the SOLVER 
implementation is still under development. This portion of the processor has been 

incorporated into an in-house system for automated vw ' lca, |°" °' ' h ® 
engineering modules. This system, the Verification Checkout Facility (VCF) 
allows interactive comparison of module results to store data as well as provides 
an intermediate checkout of the processor code. 

After validation using the VCF, the engineering modules and submodules were 
used to build a simple Technology Test Bed Engine (TTBE). The simple TT 
model has 55 states and 2 algebraic loops. This simulation of main sta 9® 
operation demonstrates steady-state and transient operation of a ,u " 
simulation, used the ROCETS general property package, and provided additoo 
information on system operation. The simple TTBE mode provides a basis for 
completion of the detailed TTBE model which will be used to verify the ROCETS 

system. 
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SYSTEM REQUIREMENT CHANGES 

The system requirements remain unchanged. 

SYSTEM ARCHITECTURE CHANGES 

character *8. The changes should increase processor sp 
search which verifies variable existence. 

COMPONENT/SUBMODULE OVERVIEW 

Modules 

Thirteen component modules have been developed lor the ROCETS system. 
R« modu^ vrere developed during Phase I of the oontraot and were used to 

mZa Cl e engine sub-system. Eight more modules w^edeve,oped p Under 

Phase II of the contract in order to model a complete rocke g 

Ses hate been slightly revised under Phase . and havener . incon»ra«ed 

into the rocket engine simulation (Simple Technology Test Bed Engine). 

The component modules fall into two general categories — state-denvative 
modules ^ non state-derivative modules. The state-derivative modu es 
perform a function that can be described by differential equations ^socia ted 
the module are "states" whose time derivatives are j calcuia' ed he 

module. Since the solver controls the value of the states. ,he mod “'“ ^ 
them as inputs and use them along with other inputs to calcu ate t 
state-derivaUves which are module outputs. Non state-derivative modules have 
“ ZZ assorted with them and are no. direcfty controlled by the solver. The 
forma for each module consists of an argument list, software identif cafton 
numter hisTory (revision dates), interface data, a description of inputs and 
outputs.' a brieT engineering description, a list of required subroutines/commons 

and finally the FORTRAN code. 

Following is a list of the component modules that have been developed or 
updated, along with a brief description. 


STATE-DERIVATIVE MODULES 


* PBRNOO 


- PREBURNER MODULE. Combined combustion and volume 
dynamics module using overall density, temperature an 
oxidizer fraction as states. 
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* MCHBOO 

* INJTOO 

* ROTROO 

* PIPEOO 

* MIXROO 

* SPLTOO 

* VOLMOO 
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MAIN COMBUSTION CHAMBER. Combustion and volume 
dynamics module using overall density, temperature and 
oxidizer fraction as states. 

MAIN FUEL INJECTOR. Perfect gas volume dynamics 
using overall density, temperature and oxidizer fraction as 
states 

- ROTOR MODULE. Simple torque balance that handles two 
supply torques, two required torques, drag torque and 
extraction torque. Rotational speed is the state. 

- PIPE WITH INERTIA. Updated pipe flow module using flow 
as a state for incompressible flow. 

- FLOW MIXER. Updated flow mixer with reverse flow 
capability using density and internal energy as states. 

FLOW SPLITTER. Updated flow splitter with reverse flow 
capability using density and internal energy as states. 

- LUMPED VOLUME. Updated volume module with reverse 
flow capability and heat transfer. Density and internal 
energy are the states. 


* PUMPOO 

* TURBOO 

* NOZLOO 

* PIPE01 

* PIPE02 
Submodules 


NON STATE-DERIVATIVE MODULES 

- PUMP MODULE. Polytropic pump with an internal iteration 
on density ratio. 

- TURBINE MODULE. Isentropic turbine with heat transfer and 
compressibility effects. 

- NOZZLE MODULE. Isentropic nozzle with normal shock. 

- PIPE WITH LOSS. Updated pipe flow module without inertia. 
Pressure loss for incompressible flow. 

FLOW THROUGH ORIFICE. Compressible flow through an 
orifice using a standard flow parameter relation. 


Twenty-one engineering submodules which are used to support the rocket 
engine component modules have been created or updated during this phase of 

the program. Basically, the sub-modules are divided into maps, utility and 
property submodules. 

Tfie map submodules, described in Table I, consist of component maps used to 
provide specific engine characteristics. 
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The utility submodules, described in Table II. consist of routines which perform 
specific functions or tasks. 


of the ROCETS 

systems. 


TABLE 1 

MAP SUBMODULES 

subroutine 

PMAP01 

description 

SSME ATD high pressure fuel pump map of both head and 
torque coefficients vs flow coefficient. 

PMAP02 

SSME ATD high pressure oxidizer (main stage) pump map 
of both head and torque coefficients vs flow coefficient. 

PMAP03 

SSME ATD high pressure oxidizer (preburner stage) pump 
map of both head and torque coefficients vs flow coefficient. 

TBMPOI 

SSME ATD high pressure oxidizer turbine map of flow 
parameter vs PR with lines of speed parameter. Also, 
efficiency over velocity ratio vs velocity ratio. 

TBMP02 

- SSME ATD high pressure fuel turbine map of flow parameter 
efficiency over velocity ratio vs velocity ratio vs PR with 
lines of speed parameter. Also, efficiency over velocity 
ratio vs velocity ratio. 

RDMP20 

- Preburner/main chamber temperature rise curves. 


TABLE II 

UTILITY SUBMODULES 

SUBROUTINE 

description 


CDNZOO 

MACH03 

# 

UNITOO 


normal shock. 

- Calculates Mach number 

_ Determines ail the Sl/English units conversion constants 
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TABLE III 

PROPERTY SUBMODULES 


SUBROUTINE 

DESCRIPTION 

COMBUSTION PROPERTIES - CREATED 

COMBOO 

- The combustion property driver routine 

PGASOO 

- The perfect gas model interface routine 

CPGMOO 

- Calculate CP and gamma for the products of combustion. 
The routine also allows for pressure ratio effects. 

CPMP01 

- CP as a function pressure & temperature for hydrogen 

RDM12 

- Molecular weight as a function of mixture ratio 

RDMP13 

- Gamma as a function of oxidizer fraction 

RDMP21 

- CP as a function of oxidizer fraction 

ZGASOO 

- Calculates compressibility (Z) effects as a function 
pressure, temperature and mixture ratio 

ZZMP01 

- Compressibility as a function of pressure and temperature 
for hydrogen 


REAL FLUID PROPERTIES - UPDATED 

PROPOO 

- Provides thermodynamic properties for the five real fluid 
properties. 


System Processors 

The configuration processor is currently under full scale development. All 
functions necessary for a single module by module checkout have been 
completed. The system can currently interpret a single module/multiple node 
input and configure the necessary FORTRAN calls. Property capability has been 
added and currently state, state derivative and additional balance capabilities are 
being coded to allow multiple module checkouts to begin. Unit definition is also 
definable through configuration input. A common structure has been defined and 
is currently in place and is described in greater detail in the architecture section. 
A linked list system is used for determining variable existence and adding new 
variables to the list during configuration. 

A* system was created called VCF (Verification Checkout Facility) which allows 
interactive comparisons to be made between a ROCETS module and a section of 
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code representing the same function in the DTM. Binary CAPA files are created 
using the DTM. This allows the engineering equations to be verified during the 
development phase. 
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Configuration Preprocessor 


Data 


Dictionary 


Configuration 

Input 


Configuration 


Preprocessor 


VCF 

Output 


Scree n 


VCF 

Input 

Screen 


Generated 


Model 


Model 


Standard Configuration Step 


off the VCF system- 




Pratt & Whitney 
FR-20282-3 


ROCETS Main System 
Configura tion Preprocessor 
First level calls from main driver routine. 
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Module 

Reading 

Routine 


Equation 

Reading 

Routine 


Property 

Reading 

Routine 


Units 

Reading 

Routine 


Externals 

Reading 

Routine 


Run 

Reading 

Routine 


States 

Reading 

Routine 
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Linked Ust 
Checking 

Routine 


Linked List 
Adding 
Routine 
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Main 

Program 

Generator 

fH 

VCFIN1 


Generate 
VCF Input 
Panel Calls 



Generate 
VCF Output 


Panel Colls 
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A example of the module input and the generated code follows. For a 
two module input, such as: • 


DEFINE SYSTEM 

MODULE: 1 USING PIPEOO , 
NODE 1 IS IN 

NODE 2 IS V0L1 ; 

MODULE: 2 USING PIPE01 
NODE 3 IS IN2 

NODE 4 IS V0L2 ; 

END SYSTEM 


The following common would have been generated by the configuration 
processor : 


COMMON / GLOBVAR 

/ 





X 

IPRPL 

p 

IUPDAT , 

M0DN1 

p 

N0D1IN 

, NGD2V0L1, 

X 

AREA1 

p 

CF1 

PTIN 

p 

PTV0L1 

, RHOIN 

X 

RH0V0L1 

p 

XLEN1 

HI 

p 

DHDT1 

, M0DN2 

X 

N0D1IN2 

p 

N0D2V0L2, 

CF2 

p 

PTIN2 

, PTV0L2 , 

X 

RH0IN2 

p 

H2 





of 

these variables 

were created 

using 

the 

standard , 

naming convention 


described in the ROCETS standards. 


The above input would also have resulted in the following calls to the 
engineering modules as listed below: 


CALL 

PIPEOO C 

IPRPL 

, IUPDAT 

, MODN1 

/ NODI IN 

p 

X 

N0D2V0L1 

, AREA1 

,CF1 

, PTIN 

, PTV0L1 

p 

X 

RHOIN 

, RHOVOL 1 

, XLEN1 

,W1 

, DHDT1 

) 

CALL 

PIPE01 ( 

IPRPL 

, IUPDAT 

f M0DN2 

/ N0D1IN2 

p 

X 

N0D2V0L2 

, CF2 

, PTIN2 

, PTV0L2* 

f RH0IN2 

p 


X W2 ) 
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The requirements for the property calls have recently been revised. 
Previously all properties were called as the following example shows: 


PROPERTY: PTV0L1 = FCRH0V0L1, UTV0L1 ) USING HYDROGEN; 

PROPERTY: PTV0L2 = FCRH0V0L2, UTV0L2 ) USING OXYGEN ; 


This would have produced the following code: 


C 

C BEGIN PROPERTY CALLS 
C 

CALL PROPOO (' HYDROGEN', 14,RH0V0L1 ,UTV0L1 ,PTV0L1 , 
X DUMYV0L1 , DUMYV0L1 , TEMP ) 

CALL PROPOO C 'OXYGEN • ,14,RH0V0L2 ,UTV0L2 ,PTV0L2 , 
X DUMYV0L2, DUMYV0L2,TEMP) 


Recently, the need for real combustion property calculations have 
altered the form of the property definition to the following: 


PROPERTY 

PROPOO 

PTV0L1 = FCRH0V0L1 , UTV0L1 ) USING HYDROGEN 
COMBOO 

<combustion definitions 

The combustion definition has not been resolved as of this printing. 
The type of calls generated would be exact like those listed above 
for PROPOO and would be similiar for COMBOO. 
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The VCF calls to the input panels from the program main for the 
above example look like: 


C 

C CHECKOUT VERIFICATION FACILITY 
C ENTRY CALL 

DO 100 1=1, 2 

IRC = 0 

DO 50 J = 1, ISIZEC I ) 

50 ARRAY! J ) = ARGLIST(I,J) 

IFCI.EQ. 2) IRC = -1 

100 CALL VCFIN C MODNM(I), MMNEM(I), ARRAY, ISIZE(I), IRC ) 


The arrays MODNM and MMENM contain the character names of the module , s 
actual subroutine name and the modules mnemonic, respectively. The 
array ARGLIST contains the argument lists as read from the data 
dictionary and ISIZE contains the actual size of each argument list. 

The variable IRC is used to verify that data was input to the simulation 
and as a single to the VCF routine that the last call has occurred. 

The VCF output panel is also invoked from the main program with the . 
following calls* 


C 

C CALL ROUTINE TO WRITE OUTPUT DATA TO VCF FILE 
C 

DO 1000 1=1, 2 

DO 500 J = 1, ISIZEC I ) 

IARRYCJ) = IPTLISTC I, J) 

500 ARRAY(J) = ARGLISTCI,J) 

1000 CALL VCFWRT ( MODNMCI), MMNEMCI), ARRAY, IARRY, ISIZE(I) ) 

C 

C INVOKE VCF OUT PANEL SYSTEM 
C 

CALL VCFOUT 

The calls to VCFWRT simply place the output from the simulation into a 
location where the VCF output routine can have access to it. The VCFOUT 
routine is the actual output panel driver. The double assignment of 
ARRAY will be removed at next update since it is redundant. 
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An example of what the VCF input screen would look like for the first 
module configured in this example: 


19 Jul 1989 = = 
Module: PIPEOO 
ROCETS 


1 

IPRPL 

2 

IUPDAT 

3 

M0DN1 

4 

NODI IN 

5 

N0D2V0L1 

6 

AREA1 

7 

CF1 

8 

PTIN 

9 

PTV0L1 

10 

RHOIN 

11 

RH0V0L1 

12 

XLEN1 

13 

Ml 

14 

DMDT1 


Verification Checkout Facility 

= = = = = = = = = = = =: = = = = = (DTM Input Map)===== =============== 3:28:37pm 

Mnemonic: 1 DTM: 048010DTM Time: 0 Page 1 

I/O DTM 


An example of what the VCF output screen format would look like for the 
first module configured in this example (without values filled in): 


Verification Checkout Facility 

19 Jul 1989 ==-=================( DTM Output Map) =================== 3:28:37pm 

Module: PIPEOO Mnemonic: 1 DTM: 048010DTM Time: 0 Page 1 

ROCETS — VAR DTM_VAR ROCETS_VAL DTM.VAL DELTA XDELTA 

1 IPRPL 

2 IUPDAT 

3 M0DN1 

4 NODI IN 

5 N0D2V0L 1 

6 AREA1 

7 CF1 

8 PTIN 

9 PTV0L1 

10 RHOIN 

11 RHOVOL1 

12 XLEN1 

13 Ml 

14 DWDT1 
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DEFINE RUN 
VCF = ON 
END RUN 

DEFINE EXTERNALS 

PTIN , HTIN , PTEX, 

CFLIN1 t HTEX , VOLVOL, 

CFLIN2; 

END EXTERNALS 

DEFINE SYSTEM 

MODULE: 1 USING PIPEOO , 

NODE 1 IS IN 

NODE 2 IS VOL 1 ; 

MODULE: 2 USING PIPE01 
NODE 3 IS IN2 

NODE <* IS V0L2 ; 

PROPERTY: PTVOL1 = F(RH0V0L1, UTVOL1 ) USING 
PROPERTY: PTV0L2 = FCRH0V0L2, UTV0L2 ) USING 

EQUATION: PT = PS * 1000/21.2 * TTIN * 

T; 

END SYSTEM 


HYDROGEN; 
OXYGEN ; 
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THIS IS THE FORTRAN THAT WOULD BE CURRENTLY PRODUCED BY THE CONFIN 
FILE LISTED ON THE PREVIOUS PAGE: 


PROGRAM MAIN 

COMMON / ARGLST / ARGIIST ( 2, 14) 

COMMON / ARGVAR / IPTLISTC 2, 14) 

COMMON / MODULE / MMNEMC 2),M0DNM( 2) 

COMMON / COMSIZ / ISIZE( 2) 

CHARACTERX8 A , ARGLIST , ARRAYC50) , MODNM , MMNEM 
INTEGER IPTLIST , IARRYC50 ) 

C 

C CHECKOUT VERIFICATION FACILITY 
C ENTRY CALL 

DO 100 1=1, 2 

IRC = 0 

DO 50 J = 1, ISIZECI) 

50 ARRAYCJ) = ARGLISTCI,J) 

IFU.EQ. 2) IRC = -1 

100 CALL VCFIN ( MODNMCI), MMNEM(I), ARRAY, ISIZE(I), IRC ) 

C 

C CALL READER ROUTINE TO BRING IN VALUES 
C 

CALL VCFREAD 
C 

C CALL MAIN SIMULATION SUBROUTINE 
C 

CALL ROCETS 
C 

C CALL ROUTINE TO WRITE OUTPUT DATA TO VCF FILE 
C 

DO 1000 1=1, 2 

DO 500 J = 1, ISIZE(I) 

IARRYCJ) = IPTLISTCI, J) 

500 ARRAYCJ) = ARGLIST(I,J) 

1000 CALL VCFWRT ( MODNMCI), MMNEMC I ) , ARRAY, IARRY, ISIZECI) ) 
C 

C INVOKE VCF OUT PANEL SYSTEM 
C 

CALL VCFOUT 
END 

SUBROUTINE ROCETS 
C 

COMMON / GLOBVAR / 


X 

IPRPL 

9 

IUPDAT , 

M0DN1 

9 

NODI IN 

9 

N0D2V0L 1 

X 

AREA1 

9 

CF1 

PTIN 

9 

PTV0L1 

9 

RHOIN 

X 

RH0V0L1 

9 

XLEN1 

HI 

9 

DHDT1 

0 

M0DN2 

X 

N0D1IN2 

9 

N0D2V0L2, 

CF2 

9 

PTIN2 

> 

PTV0L2 

X 

RH0IN2 

9 

W2 






CHARACTERX4 








X 

M0DN1 


9 NODI IN 

9 NOD2VOL1 

, N0D2V0L2 

/CF2 

X 

, PTIN2 








CALL 

UNIT DEFINITION 

ROUTINE 







C 

CALL UNITOOC 'ENGLISH ') 

C 

C BEGIN PROPERTY CALLS 
C 

CALL PROPOO C 'HYDROGEN' ,14, RH0V0L1 ,UTV0L1 ,PTV0L1 , 
* DUMYVOL 1 , DUMYVOL 1 , TEMP) 
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CALL PROPOO ('OXYGEN ',14,RH0V0L2 ,UTV0L2 
x DUMYVOL2, DUMYV0L2, TEMP ) 


, PTV0L2 


C 

c 

c 


BEGIN MODULE CALLS 


PT = PS 

CALL PIPEOO ( IPRPL 
X NOD2VOL I , AREA1 . 

x RHOIN » RHOVOL 1 


RETURN 

END 

SUBROUTINE VCFREAD 

COMMON / GLOBVAR 
X IPRPL 

X AREA1 

X RHOVOL1 , 

X N0D1IN2 , 

x RH0IN2 , 

CHARACTERX4 
X MODN1 

X , PTIN2 

NAMELIST / ROCK ■ 

X IPRPL 

X AREA1 

x RHOVOL1 

X N0DIIN2 

X RH0IN2 


X 1000/21.2 X TTIN * 
, IUPDAT 


, CF1 
, XLEN1 


, MODNl 
, PTIN 
,W1 


, NODI IN 
, PTVOL1 
, DM DTI 


IUPDAT , 

CF1 

XLEN1 

N0D2V0L2, 

M2 

, NOD1IN 


IUPDAT , 

CF1 

XLEN1 

N0D2V0L2, 

M2 


MODNl 

PTIN 

Ml 

CF2 


NODI IN 
PTVOL1 
DMDT1 
PTIN2 


MODNl 

PTIN 

Ml 

CF2 


NODI IN 
PTVOL1 
DHDT1 
PTIN2 


DO 100 1=1, 2 

100 READ (15, ROCK) 

RETURN 

END 

SUBROUTINE VCFHRT (MODNM»MMNEM, ARRAY, IARRY, ISIZE) 
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COMMON / GLOBVAR 
X IPRPL 

x AREAI 

X RHOVOL1 

X N0D1IN2 

X RH0IN2 

CHARACTERX4 
x MODNl 

X , PTIN2 


IUPDAT , 

CF1 

XLEN1 

N0D2V0L2, 

M2 

, NOD1IN 


N0D2V0L1 , 
RHOIN 
M0DN2 
PTV0L2 , 


, N0D2V0L1 , N0D2V0L2 ,CF2 


N0D2V0L1 , 
RHOIN , 
M0DN2 
PTV0L2 , 


MODNl 

t 

NODI IN , 

N0D2V0L1 

PTIN 

f 

PTVOL1 , 

RHOIN 

HI 

t 

DMDT1 

M0DN2 

CF2 

f 

PTIN2 

PTV0L2*' 

, N0D2V0L1 

, N0D2V0L2 

> * 

, CF2 


COMMON / VARTYP / TYPE( 22) 

CHARACTERX8 MO DNM , MMN EM , ARRAY ( X ) 

REALX4 VAL ( 22) 

INTEGER TYPE 


C 

INTEGER IARRY(X) 

C 

EQUIVALENCE (IPRPL ,VAL(1)) 

MRITE( 16,10) MODNM,MMNEM 
10 FORMAT! 'MODULE: ' ,A8, ' MNEMONIC: *, AS) 
DO 50 I = 1, ISIZE 
.IF (TYPE(IARRY(I)) .EQ.3) GOTO 20 
IF (TYPE(IARRYd)) .EQ.2) GOTO 30 
MRITE( 16,100) ARRAY (I),VAL(I ARRY( I ) ) 
GOTO 50 OA 


- GOTO C 5 l O b ' lZO) ARRAY(I, ' VAL ”A R RY(I), 

30 WRITE(16,140) ARRAY! I ) , VAL CIARRYf in 

50 CONTINUE «niAKKY!I)) 

100 FORMAT(3X,A8,’ = m X ,fI5.5) 

120 FORMAT! 3X, A8, ,=, ,1X,A8) 

140 F0RMAT!3X,A8,'=',l X ,I4) 

RETURN 

END 

BLOCK DATA ONE 
COMMON / GLOBNAM / A! 

COMMON / VARTYP 
COMMON / ARGLST 
COMMON ✓ ARGVAR 
COMMON / MODULE 
COMMON / COMSIZ 
CHARACTERX8 A , 

INTEGER IPTLIST 
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22 ) 

' TYPE! 22) 

' ARGLIST! 2, 14) 

' IPTLIST C 2, 14) 

' NMNEMC 2), MODNM! 2) 
' ISIZE( 2) 

ARGLIST , MMNEM , MODNM 


C 


DATA A ✓ 

' IPRPL 
'AREAl 
’RHOVOLl », 
'N0D1IN2 », 
'RH0IN2 », 


* IUPDAT 
*CF1 », 

'XLEN1 ', 

'N0D2V0L2' , 
' W2 V 


'MODNl 

’PTIN 

•W1 

•CF2 


*' 'N0D1IN •, 'N0D2V0L1*, 
'* ’PTV0L1 », 'RHOIN », 
' ’ *DWDT1 ', 'M0DN2 • , 
'» *PTIN2 •, 'PTV0L2 . . ', 


C 


C 


C 


DATA ARGLIST ✓ 

•IPRPL » IPRPL 

•M0DN2 *, 'NODI IN 
•AREAl ', 'CF2 
•PTV0L2 'PTV0L1 
•RHOVOLl », « 

’# 'DWDT1 


*' ’IUPDAT *, • IUPDAT 

’N0D1IN2 *, • N0D2V0L1 ’ , 
*' ’ CF 1 *, *PTIN2 • 

*' ' RH0IN2 «, 'RHOIN • 

'* 'XLEN1 », • , 

• . * * 


•MODNl • , 
•N0D2V0L2 ' , 
•PTIN • , 
* W2 •. 


DATA 

x 

X 

X 

X 

X 

X 


IPTLIST / 

. 1 , 

15, 

6, 

20 , 

11 , 

0, 


1, 

4, 

18, 

9, 

0, 

14, 


DATA MMNEM / 

X .j 


'2 


2 , 

16, 

7, 

21 , 

12 , 

0 / 




2 , 

5, 

19, 

10 , 

0, 


3, 

17, 

8 , 

22 , 

13, 


C 

C 


DATA MODNM ✓ 
* 'PIPEOO 


DATA ISIZE / 
* 14, 


DATA TYPE ✓ 


X 

X 

X 

X 

X 

END 


2 , 

1, 

1, 

2 , 

1, 


•PIPE01 '/ 
10/ 


2 , 

1, 

1, 

3, 

1 / 


3, 

1, 

1, 

3, 


3, 

1, 

1. 


3, 

1, 

2 , 

1, 
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Simple TTBE Simulation 

In order to complete system testing and module verification, a simple model s 
L Technology Test Bed Engine (TTBE) was completed. The a**™ 
demonstrates the ability to model a complete engine using the ROCETS 
methodology ot combining generic engineering modules with an rmpliot solven 
Figure 1 shows a schematic ot the simple TTBE along with the 42 specific stations 
in me simulation. By using generic code, only the following 13 componet modules 

are required by the simulation: 


1. INJTOO - Main Injector 

2. MCHBOO - Main Chamber 

3. MIXROO - Flow Mixer 

4. NOZLOO - Nozzle Thrust Calculations 

5. PBRNOO - Preburner 

6. PIPEOO - Incompressible flow with inertia 

7. PIPE01 - Incompresisble flow without inertia 
8 PIPE02 - Compressible flow without inertia 

9. PUMPOO - Polytropic Pump 

10. ROTROO - Rotor Torque Balance/Speed Derivative 

1 1 . SPLTOO - Flow Splitter 

12. TURBOO - Turbine 

13. VOLMOO - Volume 


The modules described above were hand processed into the simple TTBE 
simulation along with required property relationships and numerical utilities. 
There are 55 state variables and 2 algebraic loops required i" the « jMon 
as shown in Table IV. State derivatives and outputs are calculated from 
inputs and states. The simulation employs the multivariable Newton-Raphson 
soUr SMITE, to simultaneously dnve all state derivatives and algebraic loop 
errors to within a certain tolerance. SMITE operates on the matrix equation: 


AY = JAX 


Where AY is the amount that the errors need to change to be zero and AX is 
me required change in me state or algebraic loop parameter. J. the solver 
Jocabian is the matrix of partial derivatives for the states and algebraic Icraps. 
T^J matrix is evaluated by perturbing each state and 

the J matrix is inverted and the new X values are calculated until the Y variables 
are within tolerance. 

Using initial guesses from data of the Digital Transient Model P™) J at 100% 
RPL SMITE successfully drives all TTBE model state derivatives and algebras 
loop parameters to within a specific tolerance. This demonstrates the capability 
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without running rZient 00 ™ 5196 3 r ° Ckel simulaUon ,0 a steady-state point 
with small r^ P ' 9 771 simula,ion 

valve. the same ,ypa °' htovement of the fuel preburner 

valve coeffldemrat^hts^ffiw^rMotKi'm.^ theTef IbumefT' 
ret^iesh nl W3S " era,ed Un "' Ch3mber pressure (P™ 0 ^) was equaTuTthe 


fTMCHB R ,,„„ = PTMCHB„ RFL x RPL 

Xm°chb^ £ "r oxidizer fraction 

119% RPL were then run with the ^ ady T state po,nts between 60% and 

chamber pressure and LOX fraction were satis!fied n9 Th^ COeff,cients until 
of the model to use the solver ac a ^ s demonstrates the ability 

an input constraint Output of the run 'V™ 3 m ° del parameter ba sed on 

characteristic for ~ 9 

7 shows fuel and LOX preburner valve r™^ . ? Curves ' Fi 9 ure 

speeds as a function of RPL. cian,s and F ' 9ure 8 shows turbine 

C “ as a “ p ' * Fb was 

and valve data area requests calculated^ th° P COntro1 Wlth RPL re Quest input 
put through a first orSer L tl ? ? 1 °™ table ‘ The va,va ^quests were 
transients were run by inputting an RPL^en^'r'^ dynamios ' Qross throttle 

the open-loop conlro^oTSe v Ie aZ “oures^" a " d USi " 9 
run from 100% to 65% RPL rteeoi w ' p, 9 ures 9 ~13 show results of a 
109% RPL. L dGCel 3nd f,gures 14 ~ 18 ab ow a run from 65% to 
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TABLE IV - Simple TTBE Simulation States and Algebraic 

Loops 


STATE 


DESCRIPTION 

OFRFPRB 


Fuel Preburner Oxidizer Fraction 

OFRMCHB 


Main Chamber Oxidizer Fraction 

OFRMFI 


Main Fuel Injector Oxidizer Fraction 

OFROPRB 


Oxidizer Preburner Oxidizer Fraction 

RHOFMCO 


Volume FMCO Density 

RHOFPBI 


Fuel Preburner Injector Density 

RHOFPRB 


Oxidizer Preburner Injector Density 

RHOFIOS 


Volume FI OS Density 

RHOMCHB 


Main Chamber Density 

RHOMFI 


Main Fuel Injector Density 

RHOMOI 


Main Oxidizer Injector Density 

RHOOPBI 


Oxidizer Preburner Injector Density 

RHOOPRB 


Oxidizer Preburner Density 

RHOPBSF 


Preburner Fuel Splitter Density 

RHOPBSO 


Preburner Oxidizer Splitter Density 

RHOIO 


Volume 10 Density 

RH012 


Volume 12 Density 

RH015B 


Volume 15B Density 

RH08 

- 

Volume 8 Density 

SNF2 


Fuel Turbomachinery Speed 

SN02 


Oxidizer Turbomachinery Speed 

TTFPRB 


Fuel Preburner Temperature 

TTMCHB 


Main Chamber Temperature 

TTMFI 


Main Fuel Injector Temperature 

TTOPRB 


Oxidizer Preburner Temperature 

UTFMCO 


Volume FMCO Internal Energy 

UTFPBI 


Fuel Preburner Injector Internal Energy 

UTF10S 


Voume FI OS Internal Energy 

UTMOI 


Main Oxidizer Injector Internal Energy 

UTOPBI 


Oxidizer Preburner Injector Internal Energy 

UTPBSF 


Preburner Fuel Splitter Internal Energy 

UTPBSO 


Preburner Oxidizer Splitter Internal Energy 

UT10 


Volume 10 Internal Energy 

UT12 


Volume 12 Internal Energy 

UT15B 


Volume 15B Internal Energy 

UJ7S 


Volume 7S Internal Energy 

UT8 


Volume 8 Internal Energy 
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STATE 

wccv 

WFPHI 

WFPOV 

WFTI 

WFTID 

WF10 

WMFV 

WMOV 

WOPHI 

WOPOV 

W11 

W16 

W20 

W21 

W4 

W7 

W9 

ITERATION VARIABLE 


DESCRIPTION 

Coolant Control Valve Flow 

Fuel Preburner Fuel Flow 

Fuel Preburner Oxidizer Flow 

Line FTI Flow 

Line FTID Flow 

Line F10 Flow 

Main Fuel Valve Flow 

Main Oxidizer Valve Flow 

Oxidizer Preburner Fuel Flow 

Oxidizer Preburner Oxidizer Flow 

Line 11 Flow 

Line 16 Flow 

Line 20 Flow 

Line 21 Flow 

Line 4 Flow 

Line 7 Flow 

Line 9 Flow 

DESCRIPTION 


WFTRB Fuel Turbine Flow - Iterated until equal to calculated value 

WOTRB LOX Turbine flow - Iterated until equal to calculated value 
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SIMPLE TTBE SIMULATION 
FUEL PREBURNER OXIDIZER VALVE PERTUBATION 
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Figure 2 - Fuel Preburner Flow Coefficient and Flow vs Time 







SIMPLE TTBE SIMULATION 
FUEL PREBURNER OXIDIZER VALVE PERTUBATION 



Main Chamber Conditions as a Function of Time 
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TIME 

(Sec) 

FIGURE 4 - Lox Preburner Flow Coefficient and Flow vs Time 
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MAIN CHAMBER OXIDIZER FRACTION 
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Chamber Conditions as a Function of Time 
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MAIN CHAMBER OXIDIZER FRACTION 
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Rated Power Level (0-100%) 

Figure 6 - Main Chamber Conditions as a Function of RPL 
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Figure 12 - Turbine Speeds as a Funciton of Time 
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MAIN CHAMBER OXIDIZER FRACTION 
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Figure 13 - Main Chamber Conditions as a Function of Time 
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SIMPLE TTBE SIMULATION 
THROTTLE TRANSIENT ACCEL 
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as a Funciton of Time 
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Figure 17 - Turbine Speeds as a Function of Time 
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3.4.2 ENGINEERING MANUAL 

3.4.2. 1 Module Cross Reference 

3. 4.2. 2 Map Cross Reference 

3. 4. 2. 3 Map + Module Compatibility 

3. 4. 2. 4 Detailed Module Descriptions 


1. Modules 

A. Actuator 

B. Auxiliary/Special 

C. Control 

D. Heat Exchanger 

E. Injector 

1. INJTOO 

F. Line/Pipe 

1. PIPEOO 

2. PIPE01 

3. PIPE02 

G. Main Chamber 
1. MCHBOO 

H. Mixer 

1. MIXROO 

I. Nozzle 

1. NOZLOO 
j. Pre-burner 
1. PBRNOO 

K. Pump General Description 
1. PUMPOO 

L. Rotor 

1. ROTROO 

M. Secondary Flow 

N. Sensor 

O. Splitter 

1. SPLTOO 

P. Turbine 

1. TURBOO 

Q. Valve 


Page 

3.4.2. 1 

3. 4. 2. 2 

3 . 4 . 2. 3 

3. 4.2. 4 


E-1 

F-1 

F-2 

F-3 

G-1 

H-1 

1-1 

J-1 

K 

K-1 

L-1 

0-1 

P-1 
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R ' Vo,ume/ Accumulator/Tank 


1 • VOLMOO 

S. System Utility R ~ 1 

1 . SSBL03 

2. TRAP03 S-1 

3. UNITOO S-2 

2. Sub-Modules 

A. Component Sub-modules 

1. CDNZOO 

2. MACH03 2A-1 

3. FLPM02 2A-2 

Property Sub-modules 2A-3 

1 . PROPOO 

2. COMBOO 2B-1 

C. Map Sub-modules 26-2 

1- PMAPXX 

‘ 2. TBMPXX 2C-1 

D. Utility Sub-modules 2C-2 

1. SMITOI 

2. SUNB01 2D ~1 

3. PRPLOO 2D-2 

4. ITEROI 
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3. 4.2.1 


Module Cross Reference 


E-1 
F-1 
F-2 
F-3 
G-1 
H-1 
1-1 
J-1 
K-1 
L-1 
0-1 
P-1 


injtoo - 

PIPEOO 
PIPE01 
PIPE02 - 

mchboo ■ 

MIXROO 

NOZLOO 

PBBNOO 

PUMPOO 

botroo 

SPLTOO 

turboo 


Main fuel injector with inertia 

Main combustion chamber 
Simple flow mixer 
Isentropic nozzle 

Preburner without pu 9 nmress 
p'ump moduie to. polype P f ° cesS 

. Rotor 

. simple flow splitter 
- Turbine _,_ rrrank 


t 
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3. 4. 2.2 


Map Cross Reference 


1. PMAP01 

2. PMAP02 

3. PMAP03 

4. TBMP01 

5. TBMP02 


Pratt and Whitney SSME High Pressure Fuel 
Pump 

Pratt and Whitney SSME High Pressure 
Oxidizer Pump (main stage) 

Pratt and Whitney SSME High Pressure 
Oxidizer Pump (Preburner stage) 

Pratt and Whitney SSME High Pressure 
Oxidizer Turbine 

Pratt and Whitney SSME High Pressure Fuel 
Turbine 
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3. 4.2. 3 


Module + Map Compatibility 

1. PUMPOO - PMAP01, PMAP02, PMAP03 

2. TURBOO - TBMPOI . TBMP02 
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INJT00 


INPUTS: 

NODI 

NOD2 

NOD3 

NOD4 

NOD5 

NOD6 

NOD7 

NOD8 

NOD9 

NODIO 

CpF 

Cpo 

7 F 

70 ' 

7s 

71 
7e 
Q F 

Q o 

W F 

W 0 

W s 

W E 

T F 

To 

OFRp 

0FR o 


Main Fuel Injector 


| NOD2 I 
I NODI I 


W F 


| NOD5 1 
1 NOD4 | 


W 0 



Fuel Inlet Node OFR - Oxidizer Fraction 

Fuel Inlet Thermal Node T, - Injector Temperature 

Heat Flux Node - Fuel Q - Overall Density 

Oxidizer Inlet Node 

Oxidizer Inlet Thermal Node 

Heat Flux Node - Oxidizer 

Sleeve Inlet Node 

Sleeve Inlet Thermal Node 

Exit node 

Exit Thermal Node 

Specific Heat-Fuel Turbine Discharge 

Specific Heat-Oxidizer Turbine Discharge 

Specific Heat Ratio - Fuel Turbine Discharge 

Specific Heat Ratio - Oxidizer Turbine Discharge 

Specific Heat Ratio - Sleeve 

Specific Heat Ratio - Main Fuel Injector 

Specific Heat Ratio - Exit 

Heat Flux - Main Chamber Cooling - Oxidizer 

Heat Flux - Main Chamber Cooling - Fuel 


Fuel Flowrate 
Oxidizer Flowrate 
Sleeve Flowrate 
Exit Flowrate 
Fuel Temperature 
Oxidizer Temperature 
Oxidizer Fraction - Fuel Inlet 
Oxidizer Fraction - Oxidizer Inlet 
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T s - Sleeve Temperature 

T e - Exit Temperature 

Vol - Main Fuel Injector Volume 

OUTPUTS : 


d{OFR)_ _ Q X j C |j 2er Fraction Derivative 

dt 


- Density Derivative 

dt 


ip - Temperature Derivative 


CALCULATIONS: 


Win = W 5 + Wp + Wq 


Total oxygen 

flow in: w o IN = ofr f x w f + ofr 0 x w 0 


DFRIVATIVE CALCULATIONS : 

m = q x Vol 
m 

e = ~voi 

dm_{Vol) _ ( \ dVoj 

dQ _ _dt W dj 

dt {Vol) 1 

For a constant 
volume: 


dQ _ dm ( 1 \ 
~dt ~ dt\Vol) 



Win - w out 

dt 

Vol 
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ENERGY ANALYSIS: 

— = Q 0 + Qf + w f^f + W 0 h 0 + W s h s - W £ /i £ 


dU 

dt 


d(mu) du , , , dm 

— - = — (m) + («)-r- 

dt dt dt 


u = C v Tj 


— = mC v —— + C v Tj(W f + Wo + W S - We) 
dt dt 


(D 


( 2 ) 


Substitute (2) into (1) and H = CpT and T EX it = t i 

mC v — + C v Ti(W f + W 0 + W s - W E ) = Op + Oo + W F C P T F + W S C P T S - W E C P T t 
dt 

Solving for : 

< 1 Or + Oo 

^ + W 0 y o r o + W E y E Ti - (7/(Wf + W a + - W £ ))]— + — 
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BEEm 


Incompressible fluid flow 


NODI 


a 


INPUT.?- 

NODI - inlet Thermal Node 
NOD2 - Exit Thermal Node 
A - Inlet Flow Area 
CF - Flow Coefficient 
p,n - Inlet Pressure 

P°ut - Exit Pressure 

6 in - Inlet Fluid Density 

L - Length of Pipe 

OUTPUTS - 


~ Bow derivative 


£AL £|JLATI0NS- 

Definition of flow 
coefficient, OF: 


in pipe with inertia and loss 


■W- 


) (NOD2 


W 


- Flow Rate 


CFCQX J ° FJ ^ wl,ere CF - (specific gravity) 


solving for AP 


Define a : 


Qav 

Q = 
CF = 

CF 2 = 


A P = 


= l£"i 

+ Qout) 

W 

Qav 

2 

W 

j Qav ( 1 \ 

Qav 

/ Qh 2 o { A P ) 

w 2 

Qav ( 1 \ 

Qav 2 Q// 2 0 \AP J 

5 

w 2 


(1) 


a = Pav X CF 2 
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Substitute into (1): AP = — x const. (2) 

a 

Equation (2) represents the steady state relationship for W and AP. 
Lets call this friction loss AP i0SS . 

DERIVATIVE CALCULATIONS : 

Momentum analysis: 


F _ Ma 

8c 

(1) 

dV 
a = — 
dt 

(2) 

M = qAL 

(3) 

W 

W = oAV => V = 

qA 

(4) 

Substitute (2), (3) and (4) into (1): 


F _ QAL / dW \ 1 


g c \dt ) qA 


8c \dt ) 

(5) 

F = ( AP - AP l0SS ) A 

(6) 


Substitute (6) into (5) and solve for: — 

dt 


dW 

dt 


(A P - AP i oss ) 



( where AP = Pin - Pout ) 
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PIPE01 _ incompressible fluid flow in pipe with loss. 


NODI 


Q 


w 


) 


NOD2 


INPUTS: 


NODI 

Inlet Thermal Node 

NOD2 

Exit Thermal Node 

CF 

Flow Coefficient 

Pin 

Inlet Pressure 

Pout 

Exit Pressure 

Q in - 

Inlet Fluid Density 

OUTPUT: 


W - Flow Rate 


CALCULATIONS: 



Definition of flow CF __ Q x 
coefficient, CF: 

Q = — 

Qin 


where gf = — — (specific gravity) 
qh 2 o a 11 

and Q = Volumetric Flowrate (USGPM) 


[-**-( -L) 

Qin QH 2 0 \ A P } 

Solving for W: CF i _ wf. Pin 1 

p// 2 o a p 


Define a: 
Substitute into (1) 


h 2 o 


w = /ap 6in cf 2 e 

Oi = Q In X CF 2 
W = Ja A P ~x consi. 


( 1 ) 
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PIPE02 _ Compressible flow through an orifice 


NODI 


Q 


W 


3 


NOD2 


INPUTS : 

NODI 

NOD2 

K 

A 

Pin 

Pout 

R 

Tin 

7in 

OUTPUTS : 

W 


Inlet Thermal Node 
Exit Thermal Node 
Number of Head Losses 
Throat Area 
Inlet Pressure 
Exit Presure 
Gas Constant 
Inlet Temperature 
Inlet Specific Heat Ratio 


Flowrate 


CALCULATIONS : 
PR = Pin/Pout 


Submodule FLPM02 gives a value for FP from PR, K, 7, N . 
Definition of flow parameter: FP = x Tjn 

Pin x A 


Solve for W : 


w = 


FP X P/\: X A 

Jr x t in 


56 




Pratt & Whitney 
FR-20282-3 



INPUTS: 

NODI - 
NOD2 - 
NOD3 - 
NOD4 - 
NOD5 - 
NOD6 - 
NOD7 - 
NOD8 - 
OFR mf , - 
W out - 
Wfig 

W MF i - 

WmoI 

WoiG 

WoTH ~ 

VOL 

T f 

7 



STATES : 

Inlet Thermal Node OFR - Oxidizer Fraction 

Inlet Node To - Main Comb. Ch. Temperature 

Exit Node q - Overall Density 

Exit Thermal Node 

Exit Node 

Exit Thermal Node 

Main Chamber Exit Node 

Exit Thermal Node 

Oxidizer Fraction, Main Fuel Inj. 

Exit Flowrate 

Fuel Igniter Flowrate 

Main Fuel Injector Flowrate 

Main Oxid. Injector Flowrate 

Oxid. Igniter Flowrate 

Other Fuel Flowrate 

Main Chamber Volume 

Fuel Temperature 

Specific Heat Ratio (Chamber) 


diOFR) _ Q xj( jj zer f rac tj 0n derivative 

dt 
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dt Density derivative 


-—Temperature derivative 


gAL £ULATIQN.q - 

Total oxidizer 
& Fuel Flows: 

Flow-Based 
Mixture Ratio: 


m = Q x Vo l 

i 

W °IN = W MOI + W OIG + W MFI 
Wp IN = W FIG + W 0 TH + W MFI 

0Fw = w o m /Wf m 


Mass-based 

mixture ratio: 0F m = ofr/{ i _ qfr) 


Temperature Rise AT 
Where AT = /(OFR) 


?c = T F + A7 1 


L ATION.Q • 


For a 

constant volume: 


Energy analysis: 


6 = m/Vol 


dg ~~ ( Vol) _ dvp/ (m) 


at 


dt 

dt (Volf - 

dg _ Vo i ) dm 


dt 


(Vo l ) 2 


dt 


fe) 


de _ w F/n + Wom _ w> 


dt 


“out 


Vol 


dU 

dt 


~ {w f, n + w 0jN ) h • - (w out ) 


dU 

dt 


_ m du dm 

- m ~~ + u 

dt dt 


OFR mfi 

(i - ofr mf/) 


( 1 ) 

( 2 ) 
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du 

dt 


For a perfect gas 
From continuity. 


r w )h ' - (Wout) 

= I ( w Fin + W ° IN 

h = CpT - d h 55 CpdT 
O = CvT ^ du = CvdT 

m = mF + m ° 

+ VV 0|N ) - Wout 


■£>- P» 


= (Wf jn 

dt . into (3) yields', 

fort aas relationships and d( 

fluting perfect ga _ Wou( ) 

„ w CpTc ~ CvTc(WF[fj + Wo /N ^ 

. Cpr ' ' 


Substituting 


du_ 

dt 


du _ . cffi) 


dt 


dt 



OFR derivative: 


/Do __ 


/Do 


OFR - mT m 0 + mp 


Assuming 


homogeneous 


volume: (OFR) w t ex 
(1 - OFR) vV t £X 


y Vo EX 

W F EX 


d(OFR) 

dt 


_ mr 

1 

” mT 



dmt 

dt 


dmr 

dt 



/ m T 
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- (0FR> *?] 


^ = Wo - Wo = Wo w - ( OFR)Wt ex 
dt 


(D 

(2) 


dmT = ' WT EX 


( 3 ) 


dr 


Substituting (2) and (3) into (1) yields. 

d{OFR ) = _L\w 0 _ ( OFR)Wt ex - ( OFR){Wt in ~ W T EX ) ] 


dt m T L 


Simplifying: 


= _L \w 0nJ - (ofr)Wt in \ 

dt m T L /JV 
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MIXR00 - Simple Flow Mixer 


f 

Wt 

1 Q 

U 

w 2 ► 


NODI 

NOD2 

NOD3 

NOD4 


NOD5 

W 3 , 

NOD6 


INPUTS : 

NODI - Inlet 1 Thermal Node 

NOD2 - Inlet 1 Node 

NOD3 - Inlet 2 Thermal Node 

NOD4 - inlet 2 Node 

NOD5 - Exit Node 

NOD6 - Exit Thermal node 

h i - Inlet Stream 1 Enthalpy 

^2 - Inlet Stream 2 Enthalpy 

f >3 - Exit Enthalpy 

hmix - Mixer Enthalpy 

Wi - Inlet 1 Flowrate 

W 2 - Inlet 2 Flowrate 

W 3 - Exit Flowrate 

VOL - Mixer Volume 

OUTPUTS : 

— - Density Derivative 

- Internal Energy Derivative 


STATES: 

Q - Mixer Density 
u - Mixer Internal Energy 
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nFRIVATIVE QM QULATIONS: 

m = e x vol 


= W, + W 2 - W 3 

dt 


Q = 


m 


dQ_ = 
dt 


Vol 


Vo/' 


For a constant 
Volume: 


dg dm I 1 \ 
dt = dt \Vo/j 


dg 

(Wx + w 2 - Wj) 

dt 

Vo/ 


Energy: 

dU _ yy 1 hi + W 2 h 2 - W3h m i x 


dt 


dU d(mu) _ m( du\ , wf— ^ 

iu = \dt) \dt) 

Substituting (2) into (1): 


(assuming no heat transfer) 0) 

( 2 ) 


+ =W 1 h 1 +W 2 h 2 


-W 3 h m ix 


Solving for — . 


riu 

dT 


W\h\ + W2/I2 " WVimix 


/ dm\ / 1 \ 

dt )] \mj 
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NQZL00 - Isentropic Nozzle 


INPUTS : 


NODI - 
NOD2 - 
Pin 

Pout - 
Tin 

AR - 
AREA - 

7 

R 

CS - 


Inlet Thermal Node 
Exit Thermal Node 
Inlet Pressure Main Chamber 
Exit Pressure (Ambient) 

Inlet Temperature 
Area Ratio 
Throat Area 
Specific Heat Ratio 
Gas Constant 
Loss Coefficient 



OUTPUTS : 

Fg - Gross Thrust 
SI - Specific Impulse 

M e - Mach Number At Exit 

Wj - Flowrate At Throat 

CALCULATIONS : 

PR = Pin / Pout 

Submodule CDNZ00 gives values for FP and CF. 

Wt JR x Tin 

Definition of flow parameter: fp - AREA x PjN 


Solve for W t : 


w, = 


{FP) (fwHAREA) 
J(R)(T IN ) 
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r "a rea T cf\ Where CF is a thrust coefficient 

Thrust: \ Fz = Pis x area x cj 

Specific impulse: | si = f g / 
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™ N0 ~ ' Pr eburner without purge 


[nodi 1 


INPUTS- 

NODI - Oxidizer Inlet Flow Node 
NOD?' ^ uelln| et Flow Node 

Nod! I r Uel ln,6t Thermal "»*> 
Exit Flow Node 

Mixture Ratio (Light Level) 

Fuel Temperature 

Preburner Volume 

Fuel Flowrate 

Oxidizer Flowrate 

Exit Flowrate 

Specific Heat Ratio 

QUIPUIS: 

djOFR) 

dt ' 0x 'dizer fraction derivative 



OFR 

T c 

Q 


OFLT - 
T f 
VOL 
W F 
W 0 
Wqut 

7 


Oxidizer Fraction 

- Pre burner Temperature 

- Overall Density 




dt - Temperature derivative 


dg 


dt ~ Densif y derivative 


£AL £ULATlOMg : 

Flow - Based 

Mixture Ratio : of. - ^ 
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Mass - Based 

Mixture Ratio : of m = ofr / (i - ofr) 

Temperature Rise (AT): T c * = T F + AT 

Where AT = /(OFR) if preburner is lit 
or AT = const, if preburner is not lit 

DERIVATIVE CALCULATIOMS - 

6 = m/Vol 

dm_ (VOL) dVol (m) 

_ _dt ~aT 

c it ( Vol ) 2 

constant volume: ^ = ~ T7~Kr =— (t™) 

dt (Vol) 2 dt \ V0l J 

7q _ + VV 0 - I 

dt Vol 


Energy analysis: . (Wf t Wo) h . _ ^ 


dt 

dU du dm 

— = m— + w 

dt dt dt 


Substituting (2) into (1) and solving for — yields. 

dt 


du_ 

dt 


(W F + Wo)h * - (W out ) h c 


dm 1 

- »— j/„ 


For a perfect gas: ft = or => d/j = Cpdr 

u = Cv,7 =* du = CydT 


From continuity: /T? = /77p + 


rr? 0 


dm 

— = (W F + Wo) - W out 


(D 

( 2 ) 

(3) 


Substituting perfect gas relationships and — 

dt 


into equation (3) yields: 


U U r 

= l + W 0 ) C p r c - W 0 utOpT c — C v T c (W F + W 0 - W ou t) 

du_ _ djC» Tc) _ (dT c \ 
dt dt ~ Cv ^ dt ) 
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Simplifying yields: 



OFR derivative: 


Assuming 



mp 


/Dp 


OFR = 

m T 

m 0 

+ mp 


homogeneous 

volume: 

(OFR) Wt £X = 

Wo EX 


(1 

-OFR) Wr EX = 

W F E X 

d(OFR) 

dt 

(dm 

= V dt 

dmj \ 

~ m T ' d t °) 

/ m T 


= j_r 

dm 0 

dmr m o 1 



mr l 

dt 

dt mTJ 



l [ 
mr 

dm c 

dt 

L 

dmr\ 

l - (.OFR) — j 


- w 0 - Wo EX = Wo - ( OFR)Wt ex 

dt 


dmr 

HT 


= Wt in ~ Wt ex 


Substituting (2) and (3) into (1) yields. 

djOFR) _ _L[w 0 - ( OFR)Wt ex - ( OFR)(Wt in ~ VVt ex ) 
Jt m rt 


Simplifying: 


dtom. = _L [w 0 - (OFR)Wt in \ 
dt m T L 
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PI IMP00 - 


Pump Module For Polytropio Process 



INPUTS . 

NODI - Flow Node 
NOD2 - Inlet Thermal Node 
NOD3 - Exit Thermal Node 
NOD4 - Shaft Work Node 
W _ Flowrate 
P 1N - Inlet Pressure 

q |N - Inlet Density 

hlN - Inlet Enthalpy 

m _ Rotational Speed 


OUTPUTS : 


Pout - Discharge Pressure 
houT - Discharge Enthalpy 
e OUT - Discharge Density 
Tqrq - Torque Required 


.. , . n Jrn£\ until the exit pressure that results trom the 

Th,s module . era e ( ) ^ ^ pressure calculated from th e 

pSytropic equatSt ^Iterate uXpout - P-nJ- 

CALCULATIONS 1 


N„d = N X ) 

Submodule MAP gives values for 

work out _ (UQ (H ) (g£> _ 
r, = '^±h~ ( N){Torq){gc ) 


H (head) and Torque from W. e out and N. 
where gr = acceleration due to gravity 


_ (H)(gr) 
(rj) ( gc ) 


70 


Pratt & Whitney 
PR-20282-3 


houT = h, N + Ah 

A property call usinq 0 „ anW . 
at a Polytropi C process!^ ™ h ° UT 9 ' Ves a va,ue f or P, 


our. To get P OUTc lets 


look 


fo Mropic Headrise 

For a Polyiropic Process 


P. v n , c 

P • v ^ n _ q 1/n 
v = c 1/n p~1/n 


^ - headrise = J vc jp 
H = J 2 vdP 


Substitute (3) into (4) 


H = J Q 1/n p-1/n 


dP 


H »C ,/n / 2 p- 


1/n 


dP 


H = (c 1/n P 2 1_1/n _ c 1/n pw/nj pi j 
(~iTr)(c 1/n pJ-Vn_ c i/ n pi-i/nj 
Substitute (2) into (5) 

H = P 2 /n p ^~ 1/n 



-v, P^ /n pl-V nj 


(V 

( 2 ) 

(3) 

(4) 


(5) 
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For a polytropic process 

p* v n = k 



0 ) 


(2) 

(3) 


(4) 


also from 3: 



(5) 
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Pi 
Q 2 


Pj_ 
Q 1 


■■(*) 

= fl . £i 

\ n ) q 2 Q 1 

h(^) + * 

\ n J Q 1 02 
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Module implementation with e and P instead of v and P 



"h f/nI 

PoVT c = QOUt\ 

IP + QIN J 


FORTRAN implementation of p calculation 
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NODI - High Pressure Turbine Node N RO tor - Rotor Speed 

NOD2 - Low pressure Turbine node 

NOD3 - High Pressure Pump Node 

NOD4 - Low Pressure Pump Node 

T t1 - High Pressure Turbine Torque 

T X 2 - Low Pressure Turbine Torque 

T P1 - High Pressure Pump Torque 

T P2 - Low Pressure Pump Torque 

T d - Drag Torque 

T X tr - Extraction Torque 

n T1 - High Pressure Turbine Gear Ratio 

n T2 - Low Pressure Turbine Gear Ratio 

n P1 - High Pressure Pump Gear Ratio 

n P2 - Low Pressure Pump Gear Ratio 

n XTR - Extraction Gear Ratio 

l P - Polar Moment of Inertia 



(1) 


a 


dN r0 ( 0r 

dt 


(2) 


Substituting (2) 
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SPLT00 - Simple Flow Splitter 


NODI 


Wy 


NOD2 


e 

u 


NOD3 


r- VV 2 


NOD4 

1 


1 

NOD5 


vv 3 



NOD6 


INPUTS : 


STATES : 


NODI - 

Inlet Thermal Node 

NOD2 - 

Inlet Node 

NOD3 - 

Exit Node 

NOD4 - 

Exit Thermal Node 

NOD5 - 

Exit Node 

NOD6 - 

Exit Thermal Node 

h, 

Inlet Enthalpy 

h 2 

Exit Enthalpy 

h 3 

Exit Enthalpy 

hsPLT - 

Splitter Enthalpy 

VOL - 

Splitter Volume 

W , 

Inlet Flowrate 

W 2 

Exit Flowrate 

w 3 - 

Exit Flowrate 


e - Splitter Density 
u - Splitter Internal Energy 


OUTPUTS : 

^ - Density derivative 


Internal energy derivative 

dt 

CALCULATIONS : 

m = Q x vol 


— = Wt - W 2 - W 3 
dt * 

e . — 

" Vol 


mc*>o - mw 

dt Vol 2 
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For a constant 
Volume: 


dQ _ dm ( 1 
~dt ~ dt \Vol 


Tq (Wi ♦ W 2 - 
~J7 - Voi 1 


dU_ _ w n hi - W 2 h S pLT - W 3 h S PLT 

dt 


(assuming no heat transfer) CO 


dU _ d(mu) _ m(du_) + u (— ^ 

"dT ” dt \ dt ) \ dt ) 

m(du\ + u ( dm> \ = Wih, - W 2 h S PLT -W 3 h S PLT 

Substituting (2) into (1): y dt j \dt ) 


£ du , 

Solving for — • • 


= \wxht-WthsPLT ~ WihsPLT ~ u ( dt )] 
dt L 
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TURB00 - 


INPUTS : 
NODI - 
NOD2 - 
NOD3 - 
NOD4 - 
R 

7 

Z 

Pin 

Tin 

Pout 

N 

Qdot 

Qfrc 

MAP 


Turbine | nodi | 



Discharge Thermal Node 

Heat Transfer Node 

Shaft Work Node 

Inlet Gas Constant 

Inlet Specific Heat Ratio 

Inlet Compressibility Factor 

Inlet Pressure 

Inlet Temperature 

Exit Pressure 

Rotational Speed 

Heat Transfer Rate Into Turbine 

Fraction of Heat That Goes To Turbine Work 

Turbine Characteristic Map 


OUTPUTS : 

W - Turbine Flowrate 

Toro - Turbine Torque 

Tout - Discharge Temperature 

CALCULATIONS: 

PR = P in / Pout 



Submodule MAP gives values for FP, Ah' and 
Definition of Flow Parameter: FP = W J R * Tin - 


Solve for W: 


W = 


FP x P 

Jr x Tim 
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Real Ah: 


Ah = Ah' x t| 


TORQ = 


WAh 


N rad 

Exit Temperature: 


Ah = CpAT 
Cp 


= fl (^) 


Tin -Tout = 


{Ah - Qdot) 


Cp 


Tout = Tin - 


{Ah - Qdot)(,y - 


(Qdot * Qfrc) 
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VOLM00 - Lumped Volume With Energy and Continuity 

,Q 


NODI 


Win- 


NOD2 


INPUTS 


NODI - 

Inlet Thermal Node 

NOD2 - 

Inlet Node 

NOD3 - 

Exit Node 

NOD4 - 

Exit Thermal Node 

NOD5 - 

Heat Transfer Node 

Vol - 

Volume 

hin 

Inlet Enthalpy 

hout - 

Exit Enthalpy 

hvol - 

Volume Enthalpy 

Q 

Heat Flux into the Volume 

Win 

Inlet Flowrate 

Wout - 

Exit Flowrate 


r 


NOD5 


T 


Q 

U 


-Wout 


NOD4 


NOD3 


STATES : 


e 

u 


- Volume Density 

- Volume Internal Energy 


OUTPUTS : 


— - Density Derivative 

dt 


— - Internal Energy Derivative 

dt 

DFRIVATIVE CALCULATIONS: 


m = Q x vol 


dm 

~dT 


W in - Wout 


m 

6 = ~vo/ 

l dm \(Vol) _ (dVol\ (m) 

dQ _ { dt ) V dt I 

dt Vol 2 
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For a constant volume: 



dQ _ Win ~ W out 
dt ~ Vol 


Energy: 


dU_ 

dt 


- (Win) ( hin ) - ( Wout ) ( hvol ) + q 


( 1 ) 


dt/ _ d(mu) 
dt dt 



( 2 ) 


substituting (2) into (1) 



(Win) (hin) - (Wout) (hvol) +g 


solving for 


du , 

rfT' 


du 

A 


f (Win) (hin) - (Wout) (hvol) 
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SSBL03 - Steady State Balance Routine 

SSBL03 is an interface to SMIT01 which is designed to perform a steady-state 
balance on a given model. The call list for SSBL03 is defined as follows: 

ISSPNT = Print option flag 

0 = No steady-state print 

1 = Steady-state print to unit 6 
NSTATE = Number of States 
STATE = Array of State Variables 
□STATE = Array of State Variable Derivatives 

SNORM = Array of normalizing values for State Variables 

SNAME = Array of State Variable Names 

ISON = Array of on/off codes for State Variables 

0 = State not iterated 

1 = State iterated 

NVAR = Number of algebraic loops (balances) 

XVAR = Array of balance independent variables 
WAR = Array of balance dependent variables 
XNORM = Array of normalizing values for balances 
XNAME = Array of names for balances 
IXON = Array of on/off codes for balances 

0 = Balance not iterated 

1 = Balance iterated 

MXPASS = Maximum number of convergence passes 
ISAMAT = Same Jacobian Flag 

0 = Do no use previous matrix 

1 = Use previous Jacobian 

ISIG = Control flag 

0 = On first pass to set initialization 

1 = Evaluate derivatives and errors 

2 = Converged 

SSBL03 provides input to SMIT01 so that both states and algebraic loop 
(balance) parameters are iterated until their respective errors are within 
tolerance. The error terms for states are calculated as the derivative divided by 
the state and the error terms for the balance parameters are the input error terms 

(WAR). 

The states and balance parameters are normalized before going into SMIT01 by 
the input normalization values, SNORM, and XNORM. If the input normalizing 
value is zero, the initial value of the state or balance parameter will be used. 
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Values for the SMIT01 input DELTAX and DXALOW are input based on the 
normalization values. 

Specific states and balances may be turned on or off based on the input 
switching arrays ISON and IXON. SSBL03 builds the dependent and independent 
variable arrays for SMIT01 based on these inputs and the normalization 
discussed above. If the switch value is 1, the variable will be included in the 
SMIT01 array and therefore iterated, if the switch value is zero, the variable will 

not be operated on by SMIT01. 

If requested via the flag ISSPNT, SSBL03 will provide a detailed print of the 
steady-state balance. The steady-state print includes the following: 

1. The base point dependent and independent variables sent to SMIT01 

2. The dependent and independent SMIT01 variables during the generation of 
the J matrix. 

3. The dependent and independent SMIT01 variables during each convergence 
attempt along with the changes in each variable and the percent reduction in 
the sum of the squared errors. 

4. A message providing SMIT01 statistics upon completion. 
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JRAP03 - Transient Integration Routine 

sss - ssbl ° 3 and - - — ** - 


1. DT - Transient time increment 

2. TIME = Transient time 


algebraic 'toons 'Tf" "a 63 ™ 0 " r ° Uti " e WhiCh can simultaneously solve 

p r r fo r^^ nd ^^ is s ;^ to rre ^ 

e“s b be'cr;e Cin9 '° * '° SSBL ° 3 •« * 


c STATE 


SNORM 


rr„“ ZL“ “ Tm ~ 

STATE c . Mated = STATE ore>lous + (DSTATE + DSTATEp,^) X DT/2. 

SMrroi 5 te'reqSed"' der ' Va,IVe * ,h ® oalculalion °» STATE oa , cul „. a , 

soMng^onhe's^ates ^ toops a ' ,he sama «™ a as 

arrays in the same manZaTssZm*!! ** “" Struc,in9 ,he SMIT01 input 
STATES are as e STATE rather than DSTATE ~‘y 

SSBL03° r S '° r ,he al96braiC |00PS ara in the saml way as in 

ra st^U 0 a y ru^a °ni S w “" Si f ten *' 

S '! ,ail$ ' * hey ara sat W tneS"alf T R “ 

3 Z “ ~ ««— *ne value e, MXFAIL MXf”u 

TRAPoa • MXDTRD and MXFAIL ara a " internal variables to 
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allows the 
of nozzles 


A nozzle analysis requires the following input as known conditions: 

They are: PR ' AR ’ r P< " T " Ath 

|ar) ^raUo^specrtio chambeHemperature 

(To) and throat area (Ath). 

For a first test, we need to check and see if the PR is choke. The isen.rop,c tota 
to static pressure can be defined as. 


£ ■ (‘ * r -T Lmnl Y' 


(D 


At the throat, the choked PR occures at Mn - 1 

(Po\ = ( HiV* 

•Id ICHOKED \ 2 I 


ir-i 


(2) 


If PR > 



CHOKED 


then the nozzle is choked 
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„ test we need to check and see it the flow parameter at the throat 
Xr"d ne The flow parameter is de, inert as: 


FP 


_ WjTfg 



(3) 


°' ^ R °' b " " M " as 

function of PR- 

Qnivina eauation (1) tor Mn as a F(PR)- 




i K/i « a Frpm has been determined; it can be used to 
Now that a relationship tor Mn as a ( ) substitute equation (4) 

define the flow parameter as function of y and PR. Let 

into (3) 


FP = 


1 + 


©■((»■£ - ’] 
w — 

m\d(%r - ’] 



Simplifying; 



(5) 
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pressure ' he ° hoked 


FP, 


CHOKED ~ 


&rft m 


v -H . y 

2 


- 1 


]* 


"]1 

/ 2 


FP, 


CHOKED 




i/_^l 
2 


m 


y + 1 | 2 
M 


as CHOKED 

2 


Let the checked pressure term be defined 

CHOKED . j(ZZj Ljj feljSV 

FP choked = ■ CHOKED 

Let the nozzle pressure term be defined as EXIT: 


EXIT = 


••• FP 


EXIT 


= jpp^r- - ij 2 

•(&)* 


• EXIT 


Finally, 


H FPthroat > 


F F THROAT = / r / > £w • 

fp choked then the nozzle is choked. 
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UN-CHOKED nozzle 

Let's calculate the exit Mach number using subroutine Mach03 which guess 
Mn EX ij as a function of FP Exn -. 


mn EXIT ~ F(FP E xit) . 

Vary Mn EX , T until (FP Exn -) CALC - FP EXIT < tolerance 

The above balance uses subroutine ITER01 and calculates (FP EXIT ) CALC 


(FR exit) calc 


{y_ mn EXIT 

y+1 

(l + ~mn z £x»rj2(j-i) 


F P THROAT = FP EXIT * AR 



= PR 


Go calculate the thrust coefficient 
CHOKED nozzle 
If AR < 1 then 
CONVERGENT AREA 


mn E xi T = 1.0 



FP THROAT ~ FP CHOKED * AR 

Go calculate the thrust coefficient 
Else 

DIVERGENT AREA 


FP THROAT ~ FP CHOKED 


FPeXIT - FP CHOKED / AR 
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~ sarrs s — 1 a—sS" “ ~ 

Mn EX1T as a function of FP ex .t- 


Vary Mn EXtT until (FP EX it) calc 


FP ex1t < tolerance 



1st test: 

P EXIT > P AMB 


RPR = 


PR 


(*' 


jEXJT 


If RPR > 1 then 

, • *,», than ambient pressure (Pams) wM ch indicates 

The exit pressure (Pexit) ,s ^ Q calcu , ate the thrust coefticient. 

that there is no shook in the nozzle. 

2nd test Pexit 


Ps 2 = 


2ymn z £X iT _ t 


'+1 


-1 

-+1 


(after shock > Pamb) 


«-« ■ (£) ■ RFR 


If RPRS > 1, then 
coefficient. 

■ . n the nozzle Let’s guess Mn. (Mach number upstream o 
A normal shook exist m the nozzle. Le y 

the shock). 

{mnexiT ~ 21 
mn<\ = 1 + 2 

.. , D / p n across the shock to calculate. 

This allow the total to total pressure ratio (P T i t 2 ) 
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Tsl , 


P 77 

/’n 


= 1 + - — - mn 


y+\ mri\ 


i . y ~ I mn i 


y-i 


l/ts 


m/ii 


y-i 

y + 1 


r-i 


Calculating exit Mach number downstream of the shock. 


FPi 


FP 


THROAT 


EXIT ~ 


*) 


AR 


Mn Exrr = F (FP EX | T ) using subroutine Mach03 


> J EXIT l 


1 + 


7-1 )ft 

“ mn EXIT I Y 


Vary Mn! until [(Po/Ps) EX it - (PR) ( P T2 \ P T i )] < tolerance 


— | = PR 


Go calculate the thrust coefficient. 


Calculation of The Thrust Coefficient 


P 1 

*EXIT _ 1 


m 


EXIT 


The ideal throat coefficient { with P EX | T = P A mb) is defined as: 


CF, dl = y\ 


fel(rfi) 


— I y - 1 

r_1 1 1 - (\/PRyr 


The ideal thrust coefficient (no real losses, but includes expansion losses) is 
defined as: 
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CF' - FP throat mriExiT 


8 * 


(' 


r i 
2 


^lx/r) 


r 


+ AF 


p 

IJexit 

Po 



Loss real = CF'vacC 1 


Velocity coefficient is defined as: 


CV = (CF' - LOSS real) I cf idl 
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S UBMODULE FLPM02 - Flow Calculation Through A Restriction 

It is often necessary to include the effects of a flow path restriction in a simulation. 
This involves developing a relationship between flow and total pressure loss. 

FLOW CALCULATION: 

The general problem of calculating the mass flow rate through a restriction can be 
simplified considerably if the restriction is considered to be an orifice. This 
reduces the scope of the problem from developing a general theory applicable to 
all possible geometries to deriving a calculation^ procedure for only one 
geometry. The general situation is thus reduced to calculation of flow through an 
orifice and proper sizing of an orifice. 

As flow moves through an orifice it is first accelerated until it attains a maximum 
velocity at the throat. The flow is then diffused as it is decelerated to the 
downstream conditions. During this process, the total and static pressures 
change as depicted in Figure I. There is little loss in total pressure as the flow 
nears the throat, although static pressure changes considerably due to the 
increased fluid velocity. Most of the total pressure loss occurs as the fluid is 
diffused downstream of the throat. The static pressure is increasing, however, as 
a portion of the velocity head available at the throat is converted to pressure. 
Total temperature does not change throughout the process. 

Typically, the upstream and downstream total pressures and temperatures are 
known, from which flow is to be determined. Since the total pressure and 
temperature at the throat is nearly the same as upstream, flow could be 

calculated using compressible flow relations if static pressure at the throat were 
known. 

Static pressure at the throat can be calculated if it is assumed that the total 
pressure loss across the orifice is related to the velocity head at the throat. 
Numerous experiments have been conducted which justify this assumption, 
although most of the empirical data is presented in terms of discharge and flow 
coefficients. Discharge and flow coefficients are useful for control valve work but 
are generally restricted in their range of applicability. It is more convenient, 
therefore, to define a loss coefficient from which static pressure at the throat can 
be determined. 

Referring to Figure II, the loss coefficient is defined as: 

L _ Total pressure^ loss _ Pt un - Pt nnwni 
Velocity head at throat Pt up - Ps THROAT 
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Rearranging, the total to static pressure ratio at the throat is given by 


p 

s throat 


PRjL) 

1 + PR (L - 1) 



Flow 

Direction 




Figure I. Total and Static Pressure Variation in an Orifice 

where PR = P T up /Pt down and it is assumed the orifice is not choked If the 
calculated total to static pressure ratio is greater than the choke value, then 
orifice is choked, in which case, the actual total to static pressure ratio at e 

throat is equal to the choke value. 

With the total to static pressure ratio known, the flow parameter (w fRT,/PtA) can 
be determined at the throat based on compressible flow relations; 


w 


JrTt 

p t a 


FP 


JlgcY [P R —~ 0 

- ' r*7 

(y-1) PR r 



PRESSURE 
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Since total temperature remains 
assumed equal to the upstream 


constant, and total pressure at the throat 
total pressure, the mass flow rate can 


is 

be 


determined, 


w = fp 


(Ft UP ) (Aeff) 

R(T t up) 


VELOCITY HEAD AT 
THROAT 
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Total to Total Pressure Ratio (P Tup /p Tdown ) 

Flow Parameter MAP with Various Loss Coefficients 



FLOW PARAMETER WITH LOSS 
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FLOW PARAMETER MAP WITH VARIOUS LOSS COEFFICIENTS 


FLOW PARAMETER WITH LOSS 
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ROCETS REAL FLUID AND COMBUSTION PROPERTY PACKAGES 

The ROCETS property package consists of five real fluid properties and one 
"perfect gas" combustion property. Both property packages consis 
subroutine and a property/model Interface subroutine lor each ot the five fluids 
and the combustion method. The driver routine is the decision maker and I the 
routes all information to the correct property or model interface w 'hit ch 
accesses the property maps for the real fluid properties or the subroutines for 
combustion properties to look up the desired results. The program logic of the 
ROCETS real fluid properties and the combustion property package is s ow 
Figures 3.4.2.4.1 and 3A2.4.2 respectively. The five basic properties are 
para-hydrogen, oxygen, methane, nitrogen and helium. The c0 ^ bus ^ 
properttes provide a perfect gas combustion analysis using hydrogen andoxygen 
as the reactants. The driver has been set up such that it can be expande 
allow real gas combustion and Brinkly combustion at some future date. The 
combustion properties are from the digital transient model (DTM), but have been 
upgraded to meet ROCETS software standards. 

in the combustion property package, the driver receives basic inputs such as 
model type (MTYPE), property option (IOPT), the independent varies > (XVAR IT 
XVAR2, XVAR3 and XVAR4) and the map mdicies IHOLD1 and IHOLD2. 
information is then passed to the model interface which accesses the require 
information from the combustion property subroutines and maps. The c 
statement for the driver routine has the following form: 


CALL COMBOO 


(MTYPE, IOPT, XVAR1, XVAR2, XVAR3, XVAR4, OVAR1, 
OVAR2, IHOLD1, IHOLD2) 

It should be noted that MTYPE is a character's variable which defines the 
method to be used. For example, if the user wishes ^ access combus ion 
properties using a perfect gas model, they would define MTYPE = PGAS 

Once the combustion property model type has been defined, the correct model 
interface is accessed. The model interface routines have the following form. 


CALL XGASAA 


(IOPT, XVARI, XVAR2, XVAR3, XVAR4, OVAR1, OVAR2, 
IHOLD1, IHOLD2) 

The single character prefix X defines the model to be used and two character 
suffix defines the version number. The following models are available or are 
planned: P-perfect gas model (available) , R-Real gas model (not available) and 
B-Brinkly model (not available). A complete list of all options and parameters 
used in the argument list are defined In Table 3.4.2.4.4 and the units are defined 

in Table 3. 4. 2. 4. 5. 
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The perfect gas model interface defined above obtains property « n for m*.on . 
both subroutines and property maps from the Digital Trans.ent Model (DTM). 
However, the subroutines found in the DTM have been upgraded to ROCETS 
software standards and the maps which are from the DTM input file were pu in 
subroutine format. Basically, the subroutines which are called from the mode 
interface access results from tabular data which are in uni-variant or bi vana 
form. The univariant maps re set up in x-axis and y-axis coordinate sys em an 
use a linear map reader (SUNB00) to interpolate. The bi-variant maps are se 
up in an x-axis, y-axis and z-axis coordinate system and use a c ^P ond 9 
point map reader (CPMR02) to interpolate. The map indices IHOLD1 
IHOLD2 care used to track the corresponding map locations and are responsib e 
for a significant reduction in the cost of accessing the maps. 

At this stage of development, the real gas and Brinkly gas combustion models 
are not available. However, as they are implemented into the ROCETS system, 
these models will be discussed in greater detail. 

The ROCETS combustion property package has 7 property maps that are used to 
drive the perfect gas model interface. The real gas and Brinkly gas in er ace 
have not yet been implemented. The following tables define the subroutines and 

maps which are used. 

COMBUSTION PROPERTY ROUTINES 


PFRFFOT GAS REAL GAS BRINKLY GAS 

CPGM00 N0NE N0NE 

CPMP01 

RDMP13 

RDMP21 

RDMP12 

ZGAS00 

ZZMP01 

These maps are presented in figures 3.4.2.4.27 - 3.4.2.4.31. The perfect gas 
combustion property package allows the calculation of basic thermodynamic 
properties which are corrected for compressibility effects for the hot gas sys e 
components such as the preburner, turbine, injectors and the mam cham er. 

In the basic fluid property package, the driver routine receives basic property 
inputs such as fluid type (FLUID), property option (IOPT). the independent 
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variables (XVAR1 and XVAR2) and the station niumber (KSTA) . The user 
supplied information is then passed to the particular property interface which 
accesses the required information from the property maps. The call statement 
for the driver routine has the following from: 

CALL PROPOO (FLUID, IOPT, XVAR1 , XVAR2, PVAR, OVAR1 , OVAR2, KSTA) 


It should be noted that FLUID is a character* 8 variable which defines the fluid by 
name. For example, if the user whishes to access hydrogen property data, e 
would define FLUID = ‘HYDROGEN . 


Once the fluid property has been defined, the correct property interface can be 
accessed. The property interface routines for the five basic properties have the 

following form: 


CALL XXPROP (iOPT, XVAR1, XVAR2, DVAR, OVAR1 , OVAR2) 

The two character prefix XX defines the fluid to be used. The following fluids are 
available: H2 - para-hydrogen, 02 - oxygen. N2 - nitrogen, HE - helium and ME 
- methane. A complete list of all options and parameters used in the argument list 
are defined in Table 3.4.2.4.1 and units are defined in Table 3.4.2.4.2. 

Each property interface defined above obtains property data from a series of real 
fluid property maps. The call statements to the maps follow a basic format which 
uses a one letter prefix to define the fluid, a three character map variable name 
and two character version number. The map variables define how a particular 
map is set up. Let's consider an example for a hydrogen map. The fluid prefix is 
H. Hence, the call for a hydrogen map whould have the form HXYZAA. A 
hydrogen map call would have the collowing form. 


CALL HXYZAA (XX, YY, ZZ, IOPT, IHLD1, IHLD2) 

All property maps are set up in an x-axis, y-axis and z-axis coordinate system 
and uses a corresponding point map reader (CPMR02) to interpolate. e 
advantage of a corresponding point reader is that a family of discontinuous 
curves such as a property map can be adequately handled. The map reader has 
two options. Option one determines zz as a function of xx and yy. Option two 
determines xx as a function of zz and yy. The integer variables IHLD1 and IHLD2 
are map indicies which are stored internally for each call to a particular property 
interface routine by specifying KSTA. 

These indicies significantly reduces the amount of CPU time required to access 
information from the maps. 
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The ROCETS property package has 19 property data maps that are used t0 d ^ e 
the interface routines. The following table defines the maps which are available: 

REAL FLUID PROPERTY ROUTINES 

MFTHANE NITROGEN HELIU M 

MRPH01 NRHP01 ERHP01 

MTPH01 NPHT01 EPHT01 

MHPS01 NHPS01 EHPS01 

MCPT01 
MCVT01 

These maps are presented in plot form in Figures 3.4.2.4.3 - 3.4.2.4.26 and are 
restricted to the pressure and temperature ranges defined in Table 3. 4.2. 4. 3. It 
should be noted that MCPT and MCVT define Cp and Cv as a function of 
temperature with lines of constant pressure. 

This particular property package has also been optimized to minimize cost. This 
was accomplished by tracking the indices of the property maps for up to 1 00 calls 
for each fluid by storing unique integer values assigned to KSTA for each call to a 
particular fluid. A speed trial was performed which ran 100,000 calls to PROP00 
for para-hydrogen on the IBM 3083 for both an iterative and noniterative solution. 
The iterative solution required 510 seconds of real time and 1.14 CRU’s of 
execution time. The noniterative solution required 207 seconds of real time and 
0.46 CRU's of execution time. 


para-HYDROGEN OXYGEN 

HRHP01 ORHP01 

HRUP01 ORUP01 

HPUT01 OPUT01 

HHPS01 OHPS01 
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TABLE 3.4.2.4.1 


PROPERTY INTERFACE PARAMETERS OPTIONS 


IQPT 

XVAR1 

XVAR2 

1 

PRES 

TEMP 

2 

PRES 

TEMP 

4 

PRES 

ENTH 

5 

PRES 

ENTR 

6 

PRES 

ENTR 

7 

PRES 

TEMP 

8 

PRES 

TEMP 

9 

PRES 

TEMP 

10 

TEMP 

DENS 

11 

PRES 

DENS 

12 

DENS 

ENTH 

13 

PRES 

ENTH 

14 

DENS 

int.energy 

15 

PRES 

INT.ENERGY 


NOTE: Of the eleven options above, 


QVAR 

OVAR1 

OVAR2 

ENTH 

DENS 

TEMP 

ENTR 

DENS 

TEMP 

ENTR 

DENS 

TEMP 

TEMP 

DENS 

TEMP 

TEMP 

DENS 

TEMP 

DENS 

DENS 

TEMP 

CP 

GAMMA 

DENS 

TH 

CONDVISC 

DENS 

PRES 

ENTR 

ENTH 

ENTH 

TEMP 

ENTR 

PRES 

— 

* 

DENS 

— 

— 

PRES 

— 

— 

DENS 

— 


ily (9) for helium and methane are 


available. 
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TABLE 3. 4. 2. 4. 2 

PARAMETERS AND UNITS USED IN XXPROP 


PARAMETER 

UNITS 



ENGLISH 

SI 

PRESSURE 

- Lbf/IN 2 

N/m 2 

TEMPERATURE 

- DEG R 

DEG K 

DENSITY 

- LBM/IN* *3 

KG/m 3 

ENTHALPY 

- BTU/LBM 

J/Kg 

INTERNAL ENERGY 

- BTU/LBM 

J/Kg 

CP 

- BTU/IN-SEC-DEG R 

J/Kg - DEG K 

THERM. COND. 

- btu/ft-hr-r 

J/m - SEC - DEG K 

VISCOSITY 

- LBM/FT-SEC 

KG/m - SEC 

_ ENTROPY 

- BTU/LBM-DEG R 

J/Kg - DEG K 
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TABLE 3A2.4.3 


NBS PROPERTY PRESSURE AND TEMPERATURE LIMITS 


FI l IIP PROPERTY 

Para-hydrogen 

Oxygen 

Nitrogen 

Helium 

Methane 


PRFSSIJRE RANGE 
(Lbf/IN 2 ) 

0.0 - 17,404 

0.0 - 17,404 

0.0 - 20,000 

0.0 - 14,500 

1.7 - 10,000 


TfMPFRATURE RANGE 
(DEG R) 

24.8 - 720.0 

97.8 - 720.0 
113.7 - 3420.0 
5.4 - 2,700.0 
163.2 - 900.0 
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TABLE 3 . 4 . 2 . 4. 4 

COMBUSTION PROPERTY INTERFACE PARAMETER OPTIONS 


IQPT 

1 

2 

3 

4 


XVAR1 

XVAR2 

XVAR3 

PT 

TI- 

OFR 

OFR 

NA 

NA 

PT 

TT 

OFR 

RHO 

TT 

OFR 


XVAR4 

OVAR1 

OVAR2 

PR 

CP 

GAMA 

NA 

RGAS 

XMW 

NA 

Z 

NA 

NA 

PT 

RGAS 


t 
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TABLE 3. 4.2. 4.5 

PARAMETERS AND UNITS USED IN XGASAA 


PARAMETER (ARG1 


PRESSURE (PT) 

TEMPERATURE (TT) 

DENSITY (RHO) 

GAS CONSTANT (RGAS) 
SPECIFIC HEAT RATIO (GAMA) 
SPECIFIC HEAT (CP) 
COMPRESSIBILITY (Z) 
MOLECULAR WEIGHT (XMW) 
PRESSURE RATIO (PR) 
OXIDIZER FRACTION (OFR) 


UNITS 


ENGLISH 

SI 

Lbf/IN 2 

N/m 2 

°R 

°K 

Lbm/IN 3 

Kg/m 3 

Lbf IN/Lbm °R 

N/m/(Kg°K) 

D’LESS 

D’LESS 

Btu/Lbm°R) 

J/(Kg°K) 

D’LESS 

D’LESS 

Lbm/Lbm-MOLE 

Kg/Kg-MOLE 

D’LESS 

D’LESS 

D’LESS 

D’LESS- 


ft 
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TABLE 3. 4.2. 4.6 

SUBROUTINE PARAMETERS AND RANGES 


SUBROUTINE 

PARAMETER 

RANGE 

ZZMP01 

PT 

14.7 - 2000 


TT 

50 - 400 


ZZ 

0.45 - 1.55 

RDMP12 

OFR 

0 - 9999 


XMW 

2-32 

CPMP01 (1) 

PT 

50 - 5000 

TT 

100 - 2000 


PR 

1 - 4 


GAMA 

1.4 

RDMP21 

OFR 

0-1.0 


CP 

0.22 - 3.6 

RDMP13 

OFR 

0-1.0 


GAMA 

1.2 - 1.4 

NOTE: (1) 

Cp agrees with H 2 (NBS) within +/- 2% when 

turbines) is less than 

2000 PSIA and within + 


pressure is less than 5000 PSIA. 


UNITS 

Lbf/IN2 

°R 

D'LESS 

D’LESS 

Lbm/Lbm-mole 

Lbg/IN2 

°R 

D’LESS 

D’LESS 

D'LESS 

Btu/(Lbm°R) 

D'LESS 

D’LESS 

inlet pressure (for 
I- 60% when inlet 
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HYDROGEN PROPERTY MAP 
DENSITY VS PRESSURE AT CONSTANT ENTHALPY 



FIGURE 3 . 4 . 2 . 4. 4 
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TEMPERATURE - DEG R 
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PRATT & WHITNEY 
HYDROGEN PROPERTY MAP 
PRESSURE VS TEMPERATURE AT CONSTANT U 



FIGURE 3. 4. 2. 4. 5 
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PRATT & WHITNEY 
HYDROGEN PROPERTY MAP 
PRESSURE VS TEMPERATURE AT CONSTANT U 



FIGURE 3 . 4 . 2 . 4. 6 
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PRATT & WHITNEY 
HYDROGEN PROPERTY MAP 
ENTROPY VS ENTHALPY AT CONSTANT PRESSURE 



FIGURE . 3 . 4. 2 . 4. 7 
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PRATT & WHITNEY 
OXYGEN PROPERTY MAP 
DENSITY VS PRESSURE AT CONSTANT U 
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PRATT & WHITNEY 
OXYGEN PROPERTY MAP 

DENSITY VS PRESSURE AT CONSTANT ENTHALPY 
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PRATT & WHITNEY 
OXYGEN PROPERTY MAP 
PRESSURE VS TEMPERATURE AT CONSTANT U 



FIGURE 3.4.2.4.10 
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PRATT & WHITNEY 
OXYGEN PROPERTY MAP 
PRESSURE VS TEMPERATURE AT CONSTANT U 



FIGURE 3. 4. 2. 4. 1 1 
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PRATT & WHITNEY 
OXYGEN PROPERTY MAP 

ENTROPY VS ENTHALPY AT CONSTANT PRESSURE 



ENTROPY - BTU / (LBM DEG R) 

FIGURE ,3.4.2.4.12 
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PRATT & WHITNEY 
METHANE PROPERTY MAP 
DENSITY VS PRESSURE AT CONSTANT ENTHALPY 



DENSITY - LB M / FT**3 


FIGURE 3.4.2.4.13 
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PRATT & WHITNEY 
METHANE PROPERTY MAP 

ENTHALPY VS TEMPERATURE AT CONSTANT PRESSURE 



ENTHALPY - BTU / LBM 

FIGURE 3.4.2.4.14 
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PRATT & WHITNEY 
METHANE PROPERTY MAP 

ENTHALPY VS TEMPERATURE AT CONSTANT PRESSURE 



FIGURE 3.4.2.4.15 
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PRATT & WHITNEY 
METHANE PROPERTY MAP 

ENTROPY VS ENTHALPY AT CONSTANT PRESSURE 



ENTROPY - BTU / (LBM DEG R) 

FIGURE 3.4.2.4.16 
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PRATT & WHITNEY 
METHANE PROPERTY MAP 
TEMPERATURE VS CV 



FIGURE 3 . 4 . 2 . 4.17 
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PRATT & WHITNEY 
METHANE PROPERTY MAP 
TEMPERATURE VS CP 



FIGURE 3.4.2.4.18 
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PRATT & WHITNEY 
NITROGEN PROPERTY MAP 
DENSITY VS PRESSURE AT CONSTANT ENTHALPY 



FIGURE 3.4.2.4.19 
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PRATT & WHITNEY 
NITROGEN PROPERTY MAP 

PRESSURE VS TEMPERATURE AT CONSTANT ENTHALPY 



FIGURE 3.4.2.4.20 
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PRATT & WHITNEY 
NITROGEN PROPERTY MAP 

PRESSURE VS TEMPERATURE AT CONSTANT ENTHALPY 



FIGURE 3.4.2.4.21 
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PRATT & WHITNEY 
NITROGEN PROPERTY MAP 
ENTROPY VS ENTHALPY AT CONSTANT PRESSURE 



FIGURE 3.4.2.4.22 
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PRATT & WHITNEY 
HELIUM PROPERTY MAP 

DENSITY VS PRESSURE AT CONSTANT ENTHALPY 


FIGURE 3.4.2.4.23 
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PRATT & WHITNEY 
HELIUM PROPERTY MAP 

PRESSURE VS TEMPERATURE AT CONSTANT ENTHALPY 



FIGURE 3 . 4 . 2 . 4.24 
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PRATT & WHITNEY 
HELIUM PROPERTY MAP 

PRESSURE VS TEMPERATURE AT CONSTANT ENTHALPY 



PRESSURE - LBF / IN**2 
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FIGURE 3.4.2.4.25 
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PRATT & WHITNEY 
HELIUM PROPERTY MAP 

ENTROPY VS ENTHALPY AT CONSTANT PRESSURE 



LOG ENTROPY - BTU / (LBM DEG R) 

FIGURE 3.4.2.4.26 
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PRATT & WHITNEY 

COMPRESSIBLITY OF TEMPERATURE 
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Figure 3.4.2.4.27 
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PRATT & WHITNEY 

_ com BUST , on property MAP 
MOLECULAR WEIGHT VS MIXTURE RATIO 
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SPECIFIC 


PRATT & WHITNEY 

COMBUSTION PROP/fRJJc^^PRFSSURE 
HEAT VS TEMP W/LINES OF PRESbUKt 



Figure 3.4.2.4.29 
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PRATT & WHITNEY 
COMBUSTION PROPERTY MAP 
SPECIFIC HEAT VS LOX FRACTION 
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Figure 3.4.2.4.30 
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PRATT & WHITNEY 

RATin n P %?■' U r$ HR N PROPERTY MAP 
ATIO OF SPECIFIC HEATS VS LOX FRACTION 



Figure 3.4.2.4.31 
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PMAPXX - Pump Characteristic MAP Format 1 

INPUTS : 

IOPT - MAP Option 
e - MAP Density 

N - Pump Speed 

W - Flow 

SFSP - Speed Scale Factor 

SFW - Flow Scale Factor 
SFHD - Head Scale Factor 
SFTQ - Torque Scale Factor 


OUTPUTS : 


HD - Pump Head 
TORQ - Pump Torque 

DESCRIPTION OF MAP OPTION : 

(I) (IOPT = -1) : Values returned in call list are design values (100% RPL). 


(II) (IOPT=0) 


Scale factors are calculated from the input values. (HD and 
TORQ are inputs when this option is chosen). 


SFSP = N des /N 
SFHD = HD des /HD 


SFW = W DES /W 
SFTQ = TORQdes/TORQ 


(III) (IOPT = 1) 


MAP is read (curve fit) by using input values and given scale 
factors from call list. 


N map = N/SFSP 
W MAP = W/SFW 


CALCULATIONS: 


Speed (in radians): 


Tip speed is defined as: 


(MAP Speed) 
(MAP Speed) 


N r = N map /^ 


u 


p ( 2jt\ 
= 2 \ 60 y 


Nmap 


Flow coefficient (O) is 
defined as: 


<t> = 


Q_ 

AU 


(shown below) 
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- Dimensional analysis: 

By relating the speed of rotation, pump diameter and volumetric 
flow rate the following equation results: 


jjy = CONSTANT 

- (£) (£ 


We know that: 
Substituting yields: 


U a NDandA a D 2 



= CONSTANT = <t> 



Head rise coefficient can now be read from the MAP usinq flow 
coefficient and speed as inputs: 



Torque coefficient can be read from the MAP using flow coefficient as an 


NOTE: 



^ 0 L 0pti< ?? i!So T=1 ) scale factors are applied to the inputs before 
reading the MAPS as defined under Option II above. 
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PUMP HEAD: 


Dimensional analysis for pump volumetric flow rate, pressure loss, 
density and diameter yields: 


q JKp 

d 2 Q 


= CONSTANT 


= CONSTANT 

(1) 

AP 


$ 

II 

S 1 

"i j 

(2) 

UD 2 

(3) 


D z Q 

Pump head is defined as: h p = 
We know that: c 


Substituting (2) and (3) into (1) : = constant 

D*H P G r 


Hmap 


Uj_xy 

Sr 


PUMP TOROUF: 

Dimensional analysis of pump area, diameter, tip speed and density 
related to pump torque yields the following: 



For Option III (IOPT = 1) scale factors are applied to the head and torque. 


H P = Hmap * SFHD 


(pump head) 


Tp = Tmap x SFTQ 


(pump torque) 
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PMAP01 - 


Pratt & Whitney SSME High Pressure Fuel Pump 
Pump Map Format 1 
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PMAP02 - Pratt & Whitney SSME High Pressure LOX Pump 
(Main Stage) Pump Map Format 1 
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PMAP03 - Pratt & Whitney SSME High Pressure LOX Pump 
(Preburner Stage) Pump Map Format 1 
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SUBMODULE TBMPXX - Turbine Characteristic MAP Format 1 


INPUTS: 


IOPT 

N 

Tin 

P IN 

Pex 

7 

R 

SFSPD 

SFPR 

SFPSi 

SFETA 


MAP Option 

Rotational Speed 

Inlet Temperature 

Inlet Pressure 

Exit Pressure 

Specific Heat Ratio 

Gas Constant 

Speed Scale Factor 

Pressure Ratio Scale Factor 

PSI Scale Factor 

Design Point Efficiency Adder 


OUTPUTS : 


FP - Flow Parameter 

•q - Efficiency 

Ah' - Ideal Delta h 


OPTION ONE flOPT = -It 

This option will return design conditions (100% RPL) for all arguments in the call 
list. 


OPTION TWO flOPT = 01 

This option will calculate the scale factors based on the input conditions and the 
design conditions as follows: 

SFSPD = N D es/N 


SFPR 


(PRdes - D_/PRdes 
(PR - l) /PR 


SFPSI = PSI / PSI 

DES 


SFETA = (rj - r} DES ) (SFSPD ) 2 

OPTION THREE flOPT = 11 

This option reads the MAP using the input parameters and input scale factors to 
calculate the MAP parameters (actual numbers used when reading the maps). 
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CALCULATIONS ; 

Speed parameter: 


Npm = 


N 


ft X Tin 

IDEAL Ah : 

= Cp (T, n — T out) 


Ah' 

- -c P - n 

(1) 


l t in ) 

Tout 

(PoutVt- 


Tin 

V P,N ) 

(2) 

Cp 

_ yfl 



y- 1 

(3) 


Substitute (2) and (3) into (1) 



Tip Speed: 


"'T "'-T '(£)« 



Axial velocity: Axial velocity. C is obtained from the kinetic 

energy equation such that: 


energy term in the 


Ah’ = 


C 2 

2gc 


C = J2gc Ah' 


Define PSI: PSI = velocity ratio 


PSI = U/C 
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For option 3 (IOPT = 1 ) MAP parameters are calculated: 

PRMAP = PP/[PP(1 - SFPR) + SFPR] 
SNMAP = N pm x SFSPD 

PSIMAP = PSI x SFPSI 


Flow parameter can now be read from the MAP using pressure ratio and speed 
parameter as inputs: : 



Efficiency can be determined by reading t)MAP/PSIMAP as a function of PSIMAP 


tiMAP 

/PSIMAP 



EFFICIENCY: t \ = -n MAP + SFETA x (SNMAP / SNDES ) 2 
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TBMP01 - Pratt & Whitney SSME High Pressure Oxidizer Turbine 
Turbine Map Format 1. 
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TBMP02 - 


Pratt & Whitney SSME High Pressure Fuel Turbine 
Turbine Map Format 1. 
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SUBMODULE SM1T01 - Multivariable Iteration Routine 

Utility submodule SMIT01 is a multi-variable iteration routine using a modified 
Newton-Raphson technique. The matrix equation used in SMIT01 is as follows: 


ay = JAX 

Where deltaY is the amount that the errors, or dependent variables, need to 
change to be zero; and deltaX is the associated change in the independent 
variable. J, the solver Jacobian, is a matrix of partial derivative generated with the 
model. SMIT01 makes a number of passes equal to the number of iteration 
variables plus 1 through the model to generate the Jacobian. The first pass is kept 
as a base point and the each of the next n variable passes generates a column of 
the J matrix by perturbing one of the independent variables and observing the 
changes in the dependent variables from the base point. The perturbation size 
for each of the independent variables is set through the SMIT01 inpt array 
DELTAX. Each element of DELTAX is the perturbation magnitude of the 
corresponding independent variable. The J matrix is then inverted, and changes 
in the independent variable values, deltaX, are calculated as follows: 

AX = J -* AY 

The SMIT01 input array DXALOW is used to limit the change in the independent 
variables along the same direction as the deltaX vector. This is done to prevent 
excessive changes in independent variables on a single pass. This process is 
iterated until the dependent variables are within a specific tolerance. The 
maximum number of iterations allowed is input to SMIT01 as the variable 
MXPASS. The tolerance for each variable is set in the SMIT01 input array TOL. 

The arrays XVAR and YVAR contain the independent and dependent variables 
respectively. XVAR is input to SMIT01 and changed as above until each YVAR is 
within the specified tolerance. The number of variables is input as NVAR. 

To improve the efficiency of SMIT01, Broyden's algorithm for updating the 
inverse J matrix J -1 , is incorporated. Broyden’s method updates the inverse J 
matrix without evaluating or inverting a new matrix providing a large savings in 
number of total passes through the model. The matrix is updated with a secant 
method and can be performed during a convergence attempt. Broyden’s 
algorithm is applied whenever the reduction in the sum of the squared errors is 
less than the SMIT01 input SSEBRD. It may still be necessary to evaluate a new J 
matrix if SMIT01 fails to converge at an acceptable rate. If the improvement in 
the sum of the squared errors is less than SSEJAC for greater than MXFAIL 
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consecutive passes, a new J matrix is evaluated. SSEJAC and MXFAIL are 
SMIT01 inputs. 

Other SMIT01 inputs are defined as follows: 

1. IPNT = Print flag 

0 = No print 

1 = SMIT01 statistics output at completion 

2. KFILE = Logical unit number for print output 

3. ILOOP = Loop number for nesting SMIT01. Two SMIT01 loops can be 

nested; ILOOP is used to distinguish internal variable for each 
loop. 

4. ISAMAT = Flag to use previous J matrix 

0 = Do not use stored J matrix 

1 = Used stored J matrix 

5. IDTYPE = Pertubation type for J matrix calculation 

1 = Backward difference 

0 = Central difference 
+1 = Forward difference 

The flag ISIG is used to control the program flow. It is set as 0 on the first call to 
SMIT01 and is set by SMIT01 on subsequent passes as follows: 

1 = New X’s returned in an attempt to converge 

2 = Perturbation pass to create Jacobian matrix 

3 = Iteration has converged 

4 = Jacobian is singular 

5 = Maximum iterations failed to converge 

Additional parameters are output as statistics on SMIT01 operation as follows: 
PCTRDN = Percent reduction in the sum of squares error for convergence attempt 
KOUNT = Convergence attempt counter. 

KTOTAL = Total pass counter (including J matrix passes). - 
KJACOB = Number of J matrices generated 
KBROYD = Number of Boyden updates 
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SUBMODULE SUNB01 -Univariant or bi-variant sunbar-type map reader with 
option to read map in any direction and either extrapolate or return corner values. 

INPUTS 

A = Array containing map 

X = Absicssa value. Normally first independent variable 
Y = Family value. Normally second independent variable 
(Not used for univariant curve) 

T = Ordinate value. Normally dependent variable 
IOPT = OPTION CODE 

1 = X and Y are input, T is returned 

(X only for univariant) 

2 = T and Y are input, X is returned 

(T only for univariant) 

3 = X and T are input, Y is returned 

(Not applicable for univariant) 


OUTPUTS : 

IERRX = Error code describing the absicssa 

0 = Normal return 

1 = X value too small. Lowest X used. 

2 = X value too small. Extrapolated value used. 

3 = X value too large. Largest X used. 

4 = X value too large. Extrapolated value used. 

IERRY = Error code describing the family 

0 = Normal return 

1 = Y value too small. Lowest Y used. 

2 = Y value too small. Extrapolated value 'used. 

3 = Y value too large. Largest Y used. 

4 = Y value too large. Extrapolated value used. 


Method of setting up the A array 
Absicsssa = X 
Ordinate = T 
Family = Y 


A (1) = Storage space for X position indicator — initially zero 
A (2) = Storage space for Y position indicator — initially zero, 
or curve identification number 
A (3) = Number of X values defining each line 
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A (4) - Number of Y values (must be zeero for univariant curve) then, 
beginning with A(5), the map is input as: 

A (5) - 4(4+NX) = X values in ascending order 

A (5+NX) - A(4+NX+NY) = Y values in ascending order 

A (5+NX+NY) = T for X1.Y1 
A (6+NX+NY) = T for XI, Y2 

A(2*NY+NX+4) = T for X1.YNY 

A(2*NY+NX+5) = T for X2,Y1 

A(2*NY+NX+6) = T for X2, Y2 

A(3*NY+NX+4) = T for X2.YNY 

A (NX+NY * NX+N Y +4) = T for XNX.YNY 

NOTE: In the above setup, NX=number of X’S, NY=number of Y’S for bivariant 

curve. Curve for univariant, A(5) - A(4+NX) contains the X values, 
A(5+NX) - A(2*NX+4) contains the T values.. 

COMMENTS : 

1. For normal reading (IOPT=1), X is bracketed first, followed by Y bracketing. 
Since X and Y are bracketed independently, care must be exercised when 
setting up maps to ensure correct values will result. If maps are not fairly flat, 
a corresponding point map reader should be used. 

2. For reading X as a function of T and DY (IOPT=2). Y is bracketed first, 
followed by T bracketing. When corner values are returned, they look as if X 
and T were swapped. 

3. When reading Y as a function of T and X (IOPT=3), X is bracketed first, 
followed by T bracketing. This option is invalid for univariant curves. 

4. Univariant and Bivariant curves are calculated similarly but there are no 
Y-calculations. 

5. Multiple valued functions (more than 1 T for given X and Y) should not be read 
with Option 2. 


ft 


163 



